  module CollocationStatistics_Module
!------------------------------------------------------------------------------!
! Data structures and routines for CollocationStatistics.                      !
! The main data structure is Set, which contains a number of instances of File.!
!------------------------------------------------------------------------------!
  use stat, only             : Stat_Type
  use stat, only             : Init_Stat,Update_Stat,Get_Stat

  use moments, only          : Moments_Type
  use moments, only          : Init_Moments,Reset_Moments,Update_Moments,Exit_Moments

  use convert, only          : met2uv 
         
  use Compiler_Features, only: getarg_genscat, iargc_genscat

  use SWAT_Module, only      : SWAT_SelectZone

  implicit none
!------------------------------------------------------------------------------!
! Data structures.                                                             !
!------------------------------------------------------------------------------!
  type DataLine_Type
    integer                                       :: BuoyNumber
    real                                          :: BuoyLat
    real                                          :: BuoyLon
    integer                                       :: BuoyDate
    integer                                       :: BuoyTime
    real                                          :: BuoySpeed
    real                                          :: BuoyDir
    real                                          :: Buoyu
    real                                          :: buoyv
    real                                          :: ScatLat
    real                                          :: ScatLon
    integer                                       :: ScatDate
    integer                                       :: ScatTime
    real                                          :: ScatSpeed
    real                                          :: ScatDir
    real                                          :: Scatu
    real                                          :: Scatv
    real                                          :: BackSpeed
    real                                          :: BackDir
    real                                          :: Backu
    real                                          :: Backv
    logical                                       :: Common
  end type DataLine_Type

  type DataFile_Type
    integer                                       :: NrOfLines
    integer                                       :: NrOfDataLines
    integer                                       :: NrOfCommonLines
    integer                                       :: CurrentLine
    character(len=256)                            :: Name
    type(DataLine_Type), dimension(:), pointer    :: Line
  end type DataFile_Type

  type DataSet_Type
    integer                                       :: NrOfFiles
    type(DataFile_Type), dimension(:), pointer    :: File
  end type DataSet_Type

  type CollocationStatistics_Type
    integer                                       :: N
    type(Stat_Type), dimension(:), pointer        :: S
    type(Stat_Type), dimension(:), pointer        :: D
    type(Stat_Type), dimension(:), pointer        :: u
    type(Stat_Type), dimension(:), pointer        :: v
  end type CollocationStatistics_Type

  type Collocation_Type
    type(Moments_Type)                            :: u                             ! First and second moments for u
    type(Moments_Type)                            :: v                             ! First and second moments for v
  end type Collocation_Type

  type TripleCollocation_Type
    integer                                       :: N                             ! Number of files
    type(Collocation_Type), dimension(:), pointer :: TC                            ! Array of TC structs; one per file
  end type TripleCollocation_Type

  type DoubleCollocation_Type
    integer                                       :: N                             ! Number of lines
    type(Collocation_Type), dimension(:), pointer :: DC                            ! Array of DC structs; one per file
  end type DoubleCollocation_Type
!------------------------------------------------------------------------------!
! Declaration of global variables.                                             !
!------------------------------------------------------------------------------!
  integer                                         :: NrOfFiles                     ! Number of input files
  character(len=256), dimension(:), pointer       :: FileList                      ! Names of input files

  logical                                         :: DoNotClean=.false.            ! Do not remove data and multiple lines from input?
  character(len=256)                              :: Extension=' '                 ! Extension for saved input files w/o multiple entries
  character (len=4)                               :: Reference='none'              ! Reference set for triple collocation (default: no TC)
  integer                                         :: MaxIter3=0                    ! Maximum number of TC iterations (default: no TC)
  integer                                         :: MaxIter2=0                    ! Maximum number of DC iterations (default: no DC)
  logical                                         :: Common3C=.false.              ! Triple collocation only for common collocations?
  real                                            :: TCfactor=4.0                  ! Acceptance limit (default: 4-sigma)
  logical                                         :: DoBiasCorrection=.true.       ! Apply bias correction in calibration? (yes by default)
  real                                            :: r2_u=0.0, r2_v=0.0            ! Representation errors (correlation buoy-scat)
  real                                            :: csb_u=0.0, csb_v=0.0          ! Correlations scat-back (default 0)
  real                                            :: wnv_u=0.0,wnv_v=0.0           ! White noise variance in scatterometer (default 0)
  integer                                         :: BuoyIndex=1                   ! Index of buoy dataset
  integer                                         :: ScatIndex=2                   ! Index of scatterometer dataset
  integer                                         :: BackIndex=3                   ! Index of background dataset
  character(len=256)                              :: RepresentationErrorFile=' '   ! Name of file with representation errors by SCORE
  character(len=256)                              :: Results3C=' '                 ! Name of file with triple coll. results
  character(len=256)                              :: Scales3C=' '                  ! Name of file with triple coll. scalings for SCORE
  character(len=256)                              :: Errors3C=' '                  ! Name of file with triple coll. errors
  type(DataSet_Type)                              :: Set                           ! Data set
  type(CollocationStatistics_Type)                :: CommonStatistics_ScatBuoy     ! Statistics of common collocations
  type(CollocationStatistics_Type)                :: IndividualStatistics_ScatBuoy ! Statistics for each collocation file
  type(CollocationStatistics_Type)                :: NonCommonStatistics_ScatBuoy  ! Statistics of non-common statistics
  type(CollocationStatistics_Type)                :: CommonStatistics_ScatBack     ! Statistics of common collocations
  type(CollocationStatistics_Type)                :: IndividualStatistics_ScatBack ! Statistics for each collocation file
  type(CollocationStatistics_Type)                :: NonCommonStatistics_ScatBack  ! Statistics of non-common statistics
  type(TripleCollocation_Type)                    :: TripleCollocation             ! Triple collocation struct
  type(DoubleCollocation_Type)                    :: DoubleCollocation             ! Double collocation struct
  logical                                         :: StandardOnly=.false.          ! Only standard triple collocation method?

  real                                            :: MaxAngle=0.0                  ! Maximum angle between scat wind and average buoy-back
  character(len=1)                                :: Zone='E'                      ! Geographical zone selector (default: E - all Earth)
  character(len=12)                               :: Months='JFMAMJJASOND'         ! Month selector. Include all months by default
  real                                            :: MinBuoyWindSpeed=-1.0         ! Minimum buoy wind speed
  real                                            :: MaxBuoyWindSpeed=-1.0         ! Maximum buoy wind speed 

  integer                                         :: Verbosity=0
  integer, parameter                              :: DataLUN=11
  integer, parameter                              :: SaveLUN=12
  integer, parameter                              :: ResLUN=13
  integer, parameter                              :: RepLUN=14
  integer, parameter                              :: ScaleLUN=15
  integer, parameter                              :: ErrorLUN=16

  contains

!==============================================================================!
! ROUTINES FOR COMMAND LINE ARGUMENT HANDLING.                                 !
!==============================================================================!

  subroutine ReadArguments
!------------------------------------------------------------------------------!
! Read the command line arguments of CollocationStatistics.                    !
!------------------------------------------------------------------------------!
  implicit none
  integer                  :: iarg,narg,i
  character(len=256)       :: arg

  narg=iargc_genscat()
  if (narg < 1) then
    call Usage
    stop
  end if

  iarg=0
  do while (iarg < narg)
    iarg=iarg + 1
    call getarg_genscat(iarg,arg)
  
    select case (arg)
    case ('-i' , '-input')
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) NrOfFiles
      allocate(FileList(NrOfFiles))

      do i=1,NrOfFiles
        iarg=iarg + 1
        call getarg_genscat(iarg,FileList(i))
      enddo
    case ('-noclean')
      DoNotClean=.true.
    case ('-s' , '-save')
      iarg=iarg + 1
      call getarg_genscat(iarg,Extension)
    case ('-r' , '-reference')
      iarg=iarg + 1
      call getarg_genscat(iarg,Reference)
    case ('-t' , '-triple')
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) MaxIter3
    case ('-tc' , '-triplecommon')
      Common3C=.true.
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) MaxIter3
    case ('-d' , '-double')
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) MaxIter2
    case ('-f' , '-factor')
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) TCfactor
    case ('-m' , '-maxangle')
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) MaxAngle
    case ('-re' , '-representation')
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) r2_u
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) r2_v
    case ('-csb' , '-corrscatback')
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) csb_u
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) csb_v
    case ('-wnv' , 'whitenoise')
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) wnv_u
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) wnv_v
    case ('-rf' , '-reperrorfile')
      iarg=iarg + 1
      call getarg_genscat(iarg,RepresentationErrorFile)   ! Representation errors as function of s_NWP
    case ('-nobias')
      DoBiasCorrection=.false.
    case ('-z' , '-zone')
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) Zone
    case ('-months')
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) Months
    case ('wind')
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) MinBuoyWindSpeed
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) MaxBuoyWindSpeed
    case ('-res' , '-results3c')
      iarg=iarg + 1
      call getarg_genscat(iarg,Results3C)
    case ('-sc' , '-scales3c')                            ! Triple collocation scalings for use in SCORE
      iarg=iarg + 1
      call getarg_genscat(iarg,Scales3C)
    case ('-e' , '-errors3c')                             ! Triple collocation errors
      iarg=iarg + 1
      call getarg_genscat(iarg,Errors3C)
    case ('-std' , '-standardonly')
      StandardOnly=.true.
    case ('-v' , '-verbosity')
      iarg=iarg + 1
      call getarg_genscat(iarg,arg)
      read (arg,*) Verbosity
    case default
      write (*,*) 'ERROR in ReadArguments: illegal command line argument ',trim(arg)
      call Usage
      stop
   end select
  enddo
 
  if (Verbosity > 0) then
    write (*,*) ' '
    write (*,*) 'Results for routine ReadArguments:'
    write (*,*) '  File  Name'
    write (*,*) '  --------------------------------'
    write (*,'(i7,2x,a)')(i,trim(FileList(i)),i=1,NrOfFiles)
  end if

  end subroutine ReadArguments


  subroutine CheckArguments
!------------------------------------------------------------------------------!
! Checks the command line arguments.                                           !
!------------------------------------------------------------------------------!
  implicit none
  integer                  :: i
  logical                  :: FileExists

  do i=1,NrOfFiles
    inquire(file=FileList(i), Exist=FileExists)
    if (.not. FileExists) then
      write (*,*) 'ERROR in subroutine CheckArguments'
      write (*,*) 'Collocation file',i,' does not exist'
      write (*,*) 'File name is ',trim(FileList(i))
      stop
    end if
  enddo

  select case (Reference)
  case ('buoy' , 'scat' , 'back' , 'all ' , ' all' , 'none')
  case default
    write (*,*) 'ERROR in subroutine CheckArguments'
    write (*,*) 'Illegal reference parameter ',Reference
    call Usage
    stop
  end select

  if (MaxIter3 < 0) then
    write (*,*) 'ERROR in subroutine CheckArguments'
    write (*,*) 'Illegal number of triple collocation iterations',MaxIter3
    call Usage
    stop
  end if

  if (MaxIter2 < 0) then
    write (*,*) 'ERROR in subroutine CheckArguments'
    write (*,*) 'Illegal number of double collocation iterations',MaxIter2
    call Usage
    stop
  end if

  end subroutine CheckArguments


  subroutine Usage
!------------------------------------------------------------------------------!
! Prints usage.                                                                !
!------------------------------------------------------------------------------!
  implicit none

  write (*,*) ' '
  write (*,*) 'Program CollocationStatistics'
  write (*,*) ' '
  write (*,*) 'Usage: CollocationStatistics < -i | -input> <N F1 F2 ... FN> [OPTIONS]'
  write (*,*) 'with <..> : mandatory options'
  write (*,*) '     [..] : free options'
  write (*,*) '       |  : choice between alternatives'
  write (*,*) ' '
  write (*,*) 'Mandatory options:'
  write (*,*) '-i <N F1 ... FN>  : Calculate collocation statistics for N collocation'
  write (*,*) '                      files F1 to FN produced by BuoyCollocate'
  write (*,*) '-input            : Same as -i'
  write (*,*) ' '
  write (*,*) 'Free options:'
  write (*,*) '-noclean          : Do not check the input for comments and multiple'
  write (*,*) '                      entries. This option can be set when the input'
  write (*,*) '                      file has been obtained using the -s option'
  write (*,*) ' '
  write (*,*) '-s <EXT>          : Save input files without multiple entries as file'
  write (*,*) '                      FI.EXT with FI given by the -i option and EXT a'
  write (*,*) '                      character variable of length 256 at most.'
  write (*,*) '                      No files will be saved when the first character'
  write (*,*) '                      of EXT is a blank (default). Multiple entries are'
  write (*,*) '                      removed only within each file, not between files.'
  write (*,*) '-save             : Same as -s'
  write (*,*) ' '
  write (*,*) '-r <REF>          : Set REF as reference for triple collocation with'
  write (*,*) '                      REF = buoy | back | scat | none | all'
  write (*,*) '                      Default : none (no triple collocation)'
  write (*,*) '-reference        : Same as -r'
  write (*,*) ' '
  write (*,*) '-t <N3>           : Set maximum number of triple collocation iterations'
  write (*,*) '                      to N3. Default : N3 = -1 (no triple collocation)'
  write (*,*) '-triple           : Same as -t'
  write (*,*) ' '
  write (*,*) ' '
  write (*,*) '-tc <N3>          : Do triple collocation for first dataset using only'
  write (*,*) '                      collocations that are common with those of the'
  write (*,*) '                      other data sets. Set maximum number of triple'
  write (*,*) '                      collocation iterations to N3.'
  write (*,*) '                      Default : N3 = -1 (no triple collocation)'
  write (*,*) ' '
  write (*,*) '-d <N2>           : Set maximum number of double collocation iterations'
  write (*,*) '                      to N2. Default : N2 = -1 (no double collocation)'
  write (*,*) '-double           : Same as -d'
  write (*,*) ' '
  write (*,*) '-f <F>            : Omit points from collocation calculations that lie'
  write (*,*) '                      more than F times the standard deviation from the'
  write (*,*) '                      calibration curve. Default : 4.0'
  write (*,*) '-factor           : Same as -f'
  write (*,*) ' '
  write (*,*) '-m <A>            : Skip point if the difference between scatterometer'
  write (*,*) '                      wind direction and the average of buoy and'
  write (*,*) '                      background direction exceeds A. All angles in deg.'
  write (*,*) '                      This option has only effect for A > 0'
  write (*,*) '-maxangle         : Same as -m'
  write (*,*) ' '
  write (*,*) '-re <R2U R2V>     : Set representation errors R2U and R2V for the u and v'
  write (*,*) '                      components of the background (default: zero errors)'
  write (*,*) '                      Note: R2U and R2V are variances in m^2/s^2'
  write (*,*) '-representation   : Same as -re'
  write (*,*) ' '
  write (*,*) '-csb <CU CV>      : Set correlation between scatterometer and background'
  write (*,*) '                      to the values CU for u and CV for v.'
  write (*,*) '                      Default: zero correlation'
  write (*,*) '-corrscatback     : Same as -csb'
  write (*,*) ' '
  write (*,*) '-wnv <WU WV>      : Set scatterometer white noise variance to WU for u'
  write (*,*) '                      and WV for v. Default: no white noise'
  write (*,*) '-whitenoise       : Same as -wnv'
  write (*,*) ' '
  write (*,*) '-rf <F>           : Read representation errors as function of s_NWP from'
  write (*,*) '                      file F produced by SCORE'
  write (*,*) '-reperrorfile     : Same as -rf'
  write (*,*) ' '
  write (*,*) '-nobias           : Omit bias correction in triple collocation'
  write (*,*) '                      Default: include bias correction'
  write (*,*) ' '
  write (*,*) '-z <Z>            : Select geographical zone Z'
  write (*,*) '                      Z = A : all Earth (default)'
  write (*,*) '                      Z = N : Northern hemisphere'
  write (*,*) '                      Z = S : Southern hemisphere'
  write (*,*) '                      Z = T : Tropics'
  write (*,*) '                      Z = X : eXtratropics'
  write (*,*) '                      All relative to 20 degrees latitude North or South'
  write (*,*) '-zone             : Same as -z'
  write (*,*) ' '
  write (*,*) '-months <M>       : Select the months in a year to be processed (for buoy date)'
  write (*,*) '                      M should be a character string of length 12.'
  write (*,*) '                      If the character at position i equals - (minus sign)'
  write (*,*) '                      month i is skipped; otherwise it is included'
  write (*,*) '                      Default: all months included (M = JFMMAJJASOND)'
  write (*,*) '                      Example: M = JF---------D selects the winter months'
  write (*,*) ' '
  write (*,*) '-wind <A B>       : Consider only buoy wind speeds in [A,B]. A or B is neglected when < 0'
  write (*,*) ' '
  write (*,*) '-res <F>          : Write triple collocation results on file F'
  write (*,*) '-results3c        : Same as -res'
  write (*,*) ' '
  write (*,*) '-sc <F>           : Write triple collocation scalings on file F for'
  write (*,*) '                      further use as input in SCORE'
  write (*,*) '                      Note: Inactive in combination with the -rf command'
  write (*,*) '-scales3c         : Same as -sc'
  write (*,*) ' '
  write (*,*) '-std              : Apply only standard triple collocation method'
  write (*,*) '                      (Ad Stoffelen method with improved statistics)'
  write (*,*) '-standardonly     : Same as -std'
  write (*,*) ' '
  write (*,*) '-v <N>            : Set verbosity to N (default is 0: small output)'
  write (*,*) '-verbosity        : Same as -v'
  write (*,*) ' '
  write (*,*) 'Note: some triple collocation options apply only to the standard method!'

  end subroutine Usage

  
  subroutine PrintArguments
!------------------------------------------------------------------------------!
! Prints (some of) the command line arguments                                  !
!------------------------------------------------------------------------------!
  implicit none

  write (*,*) ' '
  write (*,*) ' '
  write (*,*) '============================='
  write (*,*) 'PROGRAM CollocationStatistics'
  write (*,*) '============================='
  write (*,*) ' '
  write (*,*) 'Geographical zone selected : ',Zone,'             (E: All Earth, T: Tropics, X: Extratropics,', &
              ' N: Northern hemisphere, S: Southern hemisphere)'
  write (*,*) 'Months selected            : ',Months,'  (-: skipped, any: included)'
  write (*,*) '                             JFMAMJJASOND'
  write (*,*) 'Buoy wind speed range      : ',MinBuoyWindSpeed,MaxBuoyWindSpeed

  end subroutine PrintArguments


!==============================================================================!
! ROUTINES FOR INITIALIZATION AND FINALIZATION OF DATA SET.                    !
!==============================================================================!

  subroutine InitSet
!------------------------------------------------------------------------------!
! Initializes the data set.                                                    !
!------------------------------------------------------------------------------!
  implicit none

  integer          :: NrOfLines
  integer          :: NrOfDataLines
  integer          :: i,l

  Set%NrOfFiles=NrOfFiles
  allocate(Set%File(NrOfFiles))
  do i=1,NrOfFiles
    call GetFileSize(FileList(i),NrOfLines,NrOfDataLines)  ! Find number of lines per file

    Set%File(i)%Name=FileList(i)                           ! Copy file name to set struct
    Set%File(i)%CurrentLine=0                              ! Set pointer
    Set%File(i)%NrOfLines=NrOfLines                        ! Copy number of lines per file to set struct
    Set%File(i)%NrOfDataLines=NrOfDataLines                ! Copy number of data lines per file to set struct
    Set%File(i)%NrOfCommonLines=0

    allocate(Set%File(i)%Line(NrOfDataLines))              ! Allocate data space

    do l=1,NrOfDataLines
      call InitLine(Set%File(i)%Line(l))
    enddo
  enddo

  if (Verbosity > 0) then
    write (*,*) ' '
    write (*,*) 'Results for routine InitSet:'
    write (*,*) 'Note: Multiple entries within a file are not yet removed'
    write (*,*) '       NR Of LINES'
    write (*,*) '  File Total  Data  Name'
    write (*,*) '  --------------------------------------'
    write (*,'(i7,2i6,2x,a)') (i,Set%File(i)%NrOfLines,Set%File(i)%NrOfDataLines,trim(FileList(i)),i=1,NrOfFiles)
  end if

  end subroutine InitSet


  subroutine GetFileSize(File , NrOfLines,NrOfDataLines)
!-------------------------------------------------------------------------------!
! Finds the number of data lines in a collocation file.                         !
!-------------------------------------------------------------------------------!
  implicit none

  character(len=256), intent(in)  :: File
  integer, intent(out)            :: NrOfLines
  integer, intent(out)            :: NrOfDataLines

  character(len=2)                :: Text
  integer                         :: IO

  NrOfLines=0
  NrOfDataLines=0
  open (DataLUN , file=File , status='old')
  do
    read(DataLun,*,iostat=IO) Text                         ! Read first two characters of each data line
    if (IO /= 0) then                                      ! End of file - done
      close (DataLUN, status='keep')
      return
    end if
    NrOfLines=NrOfLines + 1                                ! Increase total number of lines
    if (Text(1:1) /= '#') NrOfDataLines=NrOfDataLines + 1  ! Increase number of data lines if the line is not a comment line
  enddo

  end subroutine GetFileSize


  subroutine InitLine(Line)
!------------------------------------------------------------------------------!
! Initializes a data line.                                                     !
!------------------------------------------------------------------------------!
  implicit none
  type(DataLine_Type)  :: Line

  Line%BuoyNumber = 0
  Line%BuoyLat    = 0.0
  Line%BuoyLon    = 0.0
  Line%BuoyDate   = 0
  Line%BuoyTime   = 0
  Line%BuoySpeed  = 0.0
  Line%BuoyDir    = 0.0
  Line%Buoyu      = 0.0
  Line%Buoyv      = 0.0
  Line%ScatLat    = 0.0
  Line%ScatLon    = 0.0
  Line%ScatDate   = 0
  Line%ScatTime   = 0
  Line%ScatSpeed  = 0.0
  Line%ScatDir    = 0.0
  Line%Scatu      = 0.0
  Line%Scatv      = 0.0
  Line%BackSpeed  = 0.0
  Line%BackDir    = 0.0
  Line%Backu      = 0.0
  Line%Backv      = 0.0
  Line%Common     = .false.

  end subroutine InitLine


  subroutine ExitSet
!------------------------------------------------------------------------------!
! Deallocation.                                                                !
!------------------------------------------------------------------------------!
  implicit none
  integer  :: i

  do i=1,Set%NrOfFiles
    deallocate(Set%File(i)%Line)
  enddo
  deallocate(Set%File)

  deallocate(FileList)
 
  end subroutine ExitSet

!==============================================================================!
! ROUTINES FOR INITIALISATION AND FINALISATION OF STATISTICS STRUCTS.          !
!==============================================================================!

  subroutine InitCollocationStatistics
!------------------------------------------------------------------------------!
! Initializes structs of the CollocationStatistics_Type.                       !
!------------------------------------------------------------------------------!
  implicit none

  call InitCollocationStatisticsStruct(CommonStatistics_ScatBuoy)
  call InitCollocationStatisticsStruct(IndividualStatistics_ScatBuoy)
  call InitCollocationStatisticsStruct(NonCommonStatistics_ScatBuoy)
  call InitCollocationStatisticsStruct(CommonStatistics_ScatBack)
  call InitCollocationStatisticsStruct(IndividualStatistics_ScatBack)
  call InitCollocationStatisticsStruct(NonCommonStatistics_ScatBack)

  end subroutine InitCollocationStatistics


  subroutine InitCollocationStatisticsStruct(CollStat)
!------------------------------------------------------------------------------!
! Initializes a single struct of the CollocationStatistics_Type.               !
!------------------------------------------------------------------------------!
  implicit none
  type (CollocationStatistics_Type), intent(out) :: CollStat
  integer                                        :: i

  CollStat%N=NrOfFiles

  allocate (CollStat%S(NrOfFiles))
  allocate (CollStat%D(NrOfFiles))
  allocate (CollStat%u(NrOfFiles))
  allocate (CollStat%v(NrOfFiles))

  do i=1,NrOfFiles
    call Init_Stat(CollStat%S(i))
    call Init_Stat(CollStat%D(i))
    call Init_Stat(CollStat%u(i))
    call Init_Stat(CollStat%v(i))
  enddo

  end subroutine InitCollocationStatisticsStruct


  subroutine ExitCollocationStatistics
!------------------------------------------------------------------------------!
! Deallocates the structs of the CollocationStatistics_Type.                   !
!------------------------------------------------------------------------------!
  implicit none

  call ExitCollocationStatisticsStruct(CommonStatistics_ScatBuoy)
  call ExitCollocationStatisticsStruct(IndividualStatistics_ScatBuoy)
  call ExitCollocationStatisticsStruct(NonCommonStatistics_ScatBuoy)
  call ExitCollocationStatisticsStruct(CommonStatistics_ScatBack)
  call ExitCollocationStatisticsStruct(IndividualStatistics_ScatBack)
  call ExitCollocationStatisticsStruct(NonCommonStatistics_ScatBack)

  end subroutine ExitCollocationStatistics

 
  subroutine ExitCollocationStatisticsStruct(CollStat)
!------------------------------------------------------------------------------!
! Deallocates a single struct of the CollocationStatistics_Type.               !
!------------------------------------------------------------------------------!
  implicit none
  type (CollocationStatistics_Type), intent(inout) :: CollStat

  deallocate (CollStat%S)
  deallocate (CollStat%D)
  deallocate (CollStat%u)
  deallocate (CollStat%v)

  end subroutine ExitCollocationStatisticsStruct


  subroutine InitTripleCollocation
!------------------------------------------------------------------------------!
! Allocates and initializes the triple collocation struct.                     !
!------------------------------------------------------------------------------!
  implicit none
  integer   :: i

  TripleCollocation%N=NrOfFiles
  allocate (TripleCollocation%TC(NrOfFiles))

  do i=1,NrOfFiles
    call Init_Moments(TripleCollocation%TC(i)%u,3)
    call Init_Moments(TripleCollocation%TC(i)%v,3)
  enddo

  end subroutine InitTripleCollocation


  subroutine ExitTripleCollocation
!------------------------------------------------------------------------------!
! Deallocates the triple collocation struct.                                   !
!------------------------------------------------------------------------------!
  implicit none
  integer  :: i

  do i=1,NrOfFiles
    call Exit_Moments(TripleCollocation%TC(i)%u)
    call Exit_Moments(TripleCollocation%TC(i)%v)
  enddo

  deallocate (TripleCollocation%TC)

  end subroutine ExitTripleCollocation


  subroutine InitDoubleCollocation
!------------------------------------------------------------------------------!
! Allocates and initiates the double collocation struct.                       !
!------------------------------------------------------------------------------!
  implicit none
  integer  :: i

  DoubleCollocation%N=NrOfFiles
  allocate (DoubleCollocation%DC(NrOfFiles))

  do i=1,NrOfFiles
    call Init_Moments(DoubleCollocation%DC(i)%u,2 , ThirdOrder=.true.)
    call Init_Moments(DoubleCollocation%DC(i)%v,2 , ThirdOrder=.true.)
  enddo

  end subroutine InitDoubleCollocation


  subroutine ExitDoubleCollocation
!------------------------------------------------------------------------------!
! Deallocates the double collocation struct.                                   !
!------------------------------------------------------------------------------!
  implicit none
  integer  :: i

  do i=1,NrOfFiles
    call Exit_Moments(DoubleCollocation%DC(i)%u)
    call Exit_Moments(DoubleCollocation%DC(i)%v)
  enddo

  deallocate (DoubleCollocation%DC)

  end subroutine ExitDoubleCollocation

!==============================================================================!
! ROUTINES FOR READING INPUT DATA.                                             !
!==============================================================================!

  subroutine ReadSet
!------------------------------------------------------------------------------!
! Reads the collocation files.                                                 !
!------------------------------------------------------------------------------!
  implicit none

  type(DataLine_Type)   :: NewLine
  integer               :: i,l,N
  logical               :: Reject

  if (DoNotClean) then
!------------------------------------------------------------------------------!
! Read in all input files without checking for comment and multiple lines.     !
!------------------------------------------------------------------------------!
    do i=1,Set%NrOfFiles
      open (DataLUN , file=Set%File(i)%Name , status='old')
      Set%File(i)%CurrentLine=0
      N=0
      do l=1,Set%File(i)%NrOfLines
        call ReadLine(NewLine , Reject)               ! Read new data line from file
        call AddLine(NewLine,Set%File(i))             ! Add new data line to Set
        N=N+1                                         ! Increase counter
      enddo
      close (DataLUN , status='keep')
      Set%File(i)%NrOfDataLines=N                     ! Reset number of data lines - now without multiple entries
    enddo
  else
!------------------------------------------------------------------------------!
! Loop over all input files. Read line and check if it is already present.     !
!------------------------------------------------------------------------------!
    do i=1,Set%NrOfFiles
      open (DataLUN , file=Set%File(i)%Name , status='old')
      Set%File(i)%CurrentLine=0
      N=0
      do l=1,Set%File(i)%NrOfLines
        call ReadLine(NewLine , Reject)               ! Read new data line from file
        if (Reject) cycle                             ! New data line is a comment line - go to next line
        call CheckLine(NewLine,Set%File(i) , Reject)  ! Check new data line with those already available in Set
        if (Reject) cycle                             ! New data line is already present - go to next line
        call AddLine(NewLine,Set%File(i))             ! Add new data line to Set
        N=N+1                                         ! Increase counter
      enddo
      close (DataLUN , status='keep')
      Set%File(i)%NrOfDataLines=N                     ! Reset number of data lines - now without multiple entries
    enddo
  end if
!------------------------------------------------------------------------------!
! Print output on screen.                                                      !
!------------------------------------------------------------------------------!
  if (Verbosity > 0) then
    write (*,*) ' '
    write (*,*) 'Results for routine ReadSet:'
    write (*,*) 'Note: Multiple entries within a file are removed, so the number of data lines may have decreased'
    write (*,*) '       NR Of LINES'
    write (*,*) '  File Total  Data  Name'
    write (*,*) '  --------------------------------------'
    write (*,'(i7,2i6,2x,a)') (i,Set%File(i)%NrOfLines,Set%File(i)%NrOfDataLines,trim(FileList(i)),i=1,NrOfFiles)
  end if

  if (Verbosity > 10) then
    do i=1,Set%NrOfFiles
      write (*,*) ' '
      write (*,*) 'ReadSet: Dump of file',i,' after completion of ReadSet'
      call DumpFile(Set%File(i))
    enddo
  end if
!------------------------------------------------------------------------------!
! Save data set without multiple entries if necessary.                         !
!------------------------------------------------------------------------------!
  if (Extension(1:1) /= ' ') call SaveSet

  end subroutine ReadSet


  subroutine ReadLine(Line , Reject)
!------------------------------------------------------------------------------!
! Reads a line from a collocation file. The collocation file is assumed to be  !
! opened under logical unit number DataLUN.                                    !
!------------------------------------------------------------------------------!
  implicit none

  type (DataLine_Type), intent(out)  :: Line
  logical, intent(out)               :: Reject
  character(len=256)                 :: Text
  integer                            :: MM,m
!------------------------------------------------------------------------------!
! Read line of collocation file as characters. Reject if line starts with #.   !
!------------------------------------------------------------------------------!
  Reject=.true.
  read (DataLUN,'(a256)') Text
  if (Text(1:1) == '#') return
!------------------------------------------------------------------------------!
! Line contains data. Read data, check geographical zone and period.           !
!------------------------------------------------------------------------------!
  read (Text,*) Line%BuoyNumber,                                                                   &
                Line%BuoyLat,Line%BuoyLon,Line%BuoyDate,Line%BuoyTime,Line%BuoySpeed,Line%BuoyDir, &
                Line%ScatLat,Line%ScatLon,Line%ScatDate,Line%ScatTime,Line%ScatSpeed,Line%ScatDir, &
                                                                      Line%BackSpeed,Line%BackDir
  
  call SWAT_SelectZone(Line%BuoyLat,Zone , Reject)
  if (Reject) return

  MM=int(Line%BuoyDate/100)
  m=mod(MM,100)
  if (Months(m:m) == '-') then
    Reject=.true.
    return
  end if

  if (MinBuoyWindSpeed >= 0.0 .and. Line%BuoySpeed < MinBuoyWindSpeed) then
    Reject=.true.
    return
  end if
  if (MaxBuoyWindSpeed >= 0.0 .and. Line%BuoySpeed > MaxBuoyWindSpeed) then
    Reject=.true.
    return
  end if
!------------------------------------------------------------------------------!
! Data line within specified geographical zone. (s,d) -> (u,v).                !
!------------------------------------------------------------------------------!
  call met2uv(Line%BuoySpeed,Line%BuoyDir , Line%Buoyu,Line%Buoyv)
  call met2uv(Line%ScatSpeed,Line%ScatDir , Line%Scatu,Line%Scatv)
  call met2uv(Line%BackSpeed,Line%BackDir , Line%Backu,Line%Backv)

  end subroutine ReadLine


  subroutine CheckLine(Line,File , Reject)
!------------------------------------------------------------------------------!
! Checks if a line already is present in a file. Two lines are assumed equal   !
! if the buoy number, time, and date are the same.                             !
!------------------------------------------------------------------------------!
  implicit none

  type (DataLine_Type), intent(inout)  :: Line
  type (DataFile_Type), intent(inout)  :: File
  logical, intent(out)                 :: Reject

  integer                              :: l

  Reject=.false.
  do l=File%CurrentLine,1,-1
    if (File%Line(l)%BuoyNumber /= Line%BuoyNumber) cycle    ! Not the same buoy, next search
    if (File%Line(l)%BuoyTime   /= Line%BuoyTime)   cycle    ! Not the same buoy time, next search
    if (File%Line(l)%BuoyDate   /= Line%BuoyDate)   cycle    ! Not the same buoy date, next search

    Reject=.true.                                            ! Line already present in File, done
    return
  enddo

  end subroutine CheckLine


  subroutine AddLine(Line,File)
!------------------------------------------------------------------------------!
! Adds Line to File.                                                           !
!------------------------------------------------------------------------------!
  implicit none

  type (DataLine_Type), intent(inout)  :: Line
  type (DataFile_Type), intent(inout)  :: File

  integer                              :: l

  if (File%CurrentLine == File%NrOfDataLines) then
    write (*,*) 'ERROR in routine AddLine of module CollocationStatistics_Module:'
    write (*,*) 'Capacity exceeded in file',l,'  of dataset'
    stop
  end if

  l=File%CurrentLine + 1
  File%CurrentLine=l

  File%Line(l)%BuoyNumber = Line%BuoyNumber
  File%Line(l)%BuoyLat    = Line%BuoyLat
  File%Line(l)%BuoyLon    = Line%BuoyLon
  File%Line(l)%BuoyDate   = Line%BuoyDate
  File%Line(l)%BuoyTime   = Line%BuoyTime  
  File%Line(l)%BuoySpeed  = Line%BuoySpeed 
  File%Line(l)%BuoyDir    = Line%BuoyDir
  File%Line(l)%Buoyu      = Line%Buoyu
  File%Line(l)%Buoyv      = Line%Buoyv

  File%Line(l)%ScatLat    = Line%ScatLat   
  File%Line(l)%ScatLon    = Line%ScatLon
  File%Line(l)%ScatDate   = Line%ScatDate
  File%Line(l)%ScatTime   = Line%ScatTime
  File%Line(l)%ScatSpeed  = Line%ScatSpeed
  File%Line(l)%ScatDir    = Line%ScatDir
  File%Line(l)%Scatu      = Line%Scatu
  File%Line(l)%Scatv      = Line%Scatv

  File%Line(l)%BackSpeed  = Line%BackSpeed
  File%Line(l)%BackDir    = Line%BackDir
  File%Line(l)%Backu      = Line%Backu
  File%Line(l)%Backv      = Line%Backv
  
  end subroutine AddLine

!==============================================================================!
! ROUTINES FOR PROCESSING.                                                     !
!==============================================================================!

  subroutine FindCommonLines
!------------------------------------------------------------------------------!
! Finds the common lines in Set.                                               !
!------------------------------------------------------------------------------!
  implicit none

  integer  :: i2,l1,l2
  logical  :: Equal,Match,Search
!------------------------------------------------------------------------------!
! Check number of files in set. Print output.                                  !
!------------------------------------------------------------------------------!
  if (Set%NrOfFiles == 1) return                                               ! Nothing to do

  if (Verbosity > 2) then
    write (*,*) ' '
    write (*,*) 'CountCommonLines at start FindCommonLines:'
    call CountCommonLines
  end if
!------------------------------------------------------------------------------!
! Search for common lines (identical lines in different files).                !
!------------------------------------------------------------------------------!
  do l1=1,Set%File(1)%NrOfDataLines                                            ! Loop over all lines of File 1
    Search=.true.
    Match=.false.
    i2=1                                                                      
    do while (Search .and. i2 < Set%NrOfFiles)                                 ! Loop over File i2 with i2 from 2 to N
      i2 = i2 + 1
      Match=.false.                                                            ! Assume that File i2 does not contain Line l1 of File 1
      l2=0

      do while (.not. Match .and. l2 < Set%File(i2)%NrOfDataLines)             ! Start search in File 2
        l2=l2 + 1                                                              ! Next line of File 2
        call EqualLines(Set%File(1)%Line(l1),Set%File(i2)%Line(l2) , Equal)    ! Check if data lines are equal
        if (Equal) then                                                        ! Data lines are equal:
          Set%File(i2)%CurrentLine=l2                                          !   remember index of last data line;
          Match=.true.                                                         !   end search in File 2, continue with next file
        end if
      enddo

      if (.not. Match) Search=.false.                                          ! No match in File i2, go on with next line of File 1
    enddo

    if (Match) then                                                            ! Match found:
      Set%File(1)%Line(l1)%Common=.true.                                       !   set Line l1 in File 1 to common;
      Set%File(1)%NrOfCommonLines=Set%File(1)%NrOfCommonLines + 1              !   increase counter
      do i2=2,Set%NrOfFiles                                                    !   loop over File 2 to N:
        l2=Set%File(i2)%CurrentLine                                            !     recall index of Line in File i2 equal to Line l1,
        Set%File(i2)%Line(l2)%Common=.true.                                    !     set that Line to common,
        Set%File(i2)%NrOfCommonLines=Set%File(i2)%NrOfCommonLines + 1          !     and increase counter
      enddo
    end if
  enddo
!------------------------------------------------------------------------------!
! Print output if necessary.                                                   !
!------------------------------------------------------------------------------!
  if (Verbosity > 0) then
    write (*,*) ' '
    write (*,*) 'Results for routine FindCommonLines:'
    write (*,*) '       NR Of LINES'
    write (*,*) '  File Total  Data  Comm  Name'
    write (*,*) '  --------------------------------'
    write (*,'(i7,3i6,2x,a)') (i2,Set%File(i2)%NrOfLines,Set%File(i2)%NrOfDataLines, &
                                  Set%File(i2)%NrOfCommonLines,trim(FileList(i2)),i2=1,Set%NrOfFiles)
  end if

  if (Verbosity > 2) then
    write (*,*) ' '
    write (*,*) 'CountCommonLines at end FindCommonLines:'
    call CountCommonLines
  end if

  if (Verbosity > 11) then
    do i2=1,Set%NrOfFiles
      write (*,*) ' '
      write (*,*) 'FindCommonLines: Dump of file',i2,' after completion of FindCommonLines'
      call DumpFile(Set%File(i2))
    enddo
  end if

  end subroutine FindCommonLines


  subroutine EqualLines(Line1,Line2 , Equal)
!------------------------------------------------------------------------------!
! Check if two lines are equal. The two lines are assumed equal if the         !
! following attributes are the same: BuoyNumber, BuoyDate, BuoyTime,           !
! ScatDate, and ScatTime.                                                      !
!------------------------------------------------------------------------------!
  implicit none

  type(DataLine_Type), intent(in)  :: Line1,Line2
  logical, intent(out)             :: Equal

  Equal=.false.
  if (Line1%BuoyNumber /= Line2%BuoyNumber) return
  if (Line1%BuoyDate   /= Line2%BuoyDate)   return
  if (Line1%BuoyTime   /= Line2%BuoyTime)   return
  
  if (Line1%ScatDate   /= Line2%ScatDate)   return
  if (Line1%ScatTime   /= Line2%ScatTime)   return
  Equal=.true.

  end subroutine EqualLines


  subroutine CountCommonLines
!------------------------------------------------------------------------------!
! Counts the number of common and non-common lines in a file.                  !
!------------------------------------------------------------------------------!
  implicit none

  integer    :: i,l
  integer    :: Total,Common,NonCommon

  write (*,*) 'Results from CountCommonLines:'
  write (*,*) '  File     Total    Common NonCommon'
  write (*,*) '  ----------------------------------'

  do i=1,Set%NrOfFiles
    Total=0
    Common=0
    NonCommon=0
    do l=1,Set%File(i)%NrOfDataLines
      Total=Total + 1
      if (Set%File(i)%Line(l)%Common) then
        Common=Common + 1
      else
        NonCommon=NonCommon + 1
      end if
    enddo

    write (*,'(i7,3i10)') i,Total,Common,NonCommon
  enddo

  end subroutine CountCommonLines
  

  subroutine CalcCollocationStatistics
!------------------------------------------------------------------------------!
! Calculates the collocation statistics.                                       !
!------------------------------------------------------------------------------!
  implicit none

  real     :: ObsS,ObsD,Obsu,Obsv
  real     :: SelS,SelD,Selu,Selv
  real     :: dS,dD,du,dv
  integer  :: i,l

  do i=1,Set%NrOfFiles                                  ! Loop over all collocation files
    do l=1,Set%File(i)%NrOfDataLines                    ! Loop over all collocations in file i
      SelS=Set%File(i)%Line(l)%ScatSpeed                ! Get scatterometer speed and direction
      SelD=Set%File(i)%Line(l)%ScatDir
      Selu=Set%File(i)%Line(l)%Scatu                    ! Get scatterometer u and v
      Selv=Set%File(i)%Line(l)%Scatv
!------------------------------------------------------------------------------!
!     Buoy minus scatterometer.                                                !
!------------------------------------------------------------------------------!
      ObsS=Set%File(i)%Line(l)%BuoySpeed                ! Get buoy speed and direction
      ObsD=Set%File(i)%Line(l)%BuoyDir
      Obsu=Set%File(i)%Line(l)%Buoyu                    ! Get buoy u and v
      Obsv=Set%File(i)%Line(l)%Buoyv

      dS=ObsS - SelS                                    ! Calculate differences
      dD=ObsD - SelD
      du=Obsu - Selu
      dv=Obsv - Selv

      if (dD < -180.0) dD=dD + 360.0                    ! Make sure that difference in direction lies in (-180,180)
      if (dD >  180.0) dD=dD - 360.0

      call Update_Stat(IndividualStatistics_ScatBuoy%S(i) , dS)  ! Update the statistics for all collocations in a file
      call Update_Stat(IndividualStatistics_ScatBuoy%D(i) , dD)
      call Update_Stat(IndividualStatistics_ScatBuoy%u(i) , du) 
      call Update_Stat(IndividualStatistics_ScatBuoy%v(i) , dv)

      if (Set%File(i)%Line(l)%Common) then
        call Update_Stat(CommonStatistics_ScatBuoy%S(i) , dS)    ! Update the statistics for the common collocations
        call Update_Stat(CommonStatistics_ScatBuoy%D(i) , dD)
        call Update_Stat(CommonStatistics_ScatBuoy%u(i) , du)
        call Update_Stat(CommonStatistics_ScatBuoy%v(i) , dv)
      else
        call Update_Stat(NonCommonStatistics_ScatBuoy%S(i) , dS) ! Update the statistics for the non-common collocations
        call Update_Stat(NonCommonStatistics_ScatBuoy%D(i) , dD)
        call Update_Stat(NonCommonStatistics_ScatBuoy%u(i) , du)
        call Update_Stat(NonCommonStatistics_ScatBuoy%v(i) , dv)
      end if
!------------------------------------------------------------------------------!
!     Background minus scatterometer.                                          !
!------------------------------------------------------------------------------!
      ObsS=Set%File(i)%Line(l)%BackSpeed                ! Get model speed and direction
      ObsD=Set%File(i)%Line(l)%BackDir
      Obsu=Set%File(i)%Line(l)%Backu                    ! Get model u and v
      Obsv=Set%File(i)%Line(l)%Backv

      dS=ObsS - SelS                                    ! Calculate differences
      dD=ObsD - SelD
      du=Obsu - Selu
      dv=Obsv - Selv

      if (dD < -180.0) dD=dD + 360.0                    ! Make sure that difference in direction lies in (-180,180)
      if (dD >  180.0) dD=dD - 360.0

      call Update_Stat(IndividualStatistics_ScatBack%S(i) , dS)  ! Update the statistics for all collocations in a file
      call Update_Stat(IndividualStatistics_ScatBack%D(i) , dD)
      call Update_Stat(IndividualStatistics_ScatBack%u(i) , du)
      call Update_Stat(IndividualStatistics_ScatBack%v(i) , dv)

      if (Set%File(i)%Line(l)%Common) then
        call Update_Stat(CommonStatistics_ScatBack%S(i) , dS)    ! Update the statistics for the common collocations
        call Update_Stat(CommonStatistics_ScatBack%D(i) , dD)
        call Update_Stat(CommonStatistics_ScatBack%u(i) , du)
        call Update_Stat(CommonStatistics_ScatBack%v(i) , dv)
      else
        call Update_Stat(NonCommonStatistics_ScatBack%S(i) , dS) ! Update the statistics for the non-common collocations
        call Update_Stat(NonCommonStatistics_ScatBack%D(i) , dD)
        call Update_Stat(NonCommonStatistics_ScatBack%u(i) , du)
        call Update_Stat(NonCommonStatistics_ScatBack%v(i) , dv)
      end if
    enddo
  enddo

  end subroutine CalcCollocationStatistics


  subroutine CalcTripleCollocation
!------------------------------------------------------------------------------!
! Calculate triple collocation.                                                !
!------------------------------------------------------------------------------!
  implicit none
  integer                        :: i
  real, dimension(3)             :: eps_u,eps_v
  real, dimension(3)             :: Scale_u,Scale_v
  character(len=256)             :: In1,In2
  character(len=3)               :: Par1,Par2
  real                           :: Resol
  integer                        :: N,Nfour
  real                           :: Scale,r2_t,r2_l
  integer                        :: Nset
!------------------------------------------------------------------------------!
! Check number of iterations. Open file with results. Set background index.    !
!------------------------------------------------------------------------------!
  if (MaxIter3 == 0) return

  if (Results3C(1:1) /= ' ') open (ResLUN , file=Results3C , status='unknown')

  select case (Reference)
  case ('buoy')
    BuoyIndex=1
    ScatIndex=2
    BackIndex=3
  case ('scat')
    ScatIndex=1
    BackIndex=2
    BuoyIndex=3
  case ('back')
    BackIndex=1
    BuoyIndex=2
    ScatIndex=3
  case default
    write (*,*) 'ERROR in CalcTripleCollocation: unknown reference ',Reference
    stop
  end select

  write (*,*) 'CalcTripleCollocation entered with reference ',Reference
!------------------------------------------------------------------------------!
! No file with representation errors as function of s_NWP given:               !
! - do triple collocation for first set only if Common3C == .true. or          !
! - compare methods for all data sets unless StandardOnly == true.             !
!------------------------------------------------------------------------------!
  if (RepresentationErrorFile(1:1) == ' ') then
    Nset=Set%NrOfFiles
    if (Common3C) Nset=1
    do i=1,Nset
      if (.not. StandardOnly) call CalcTCforFile_Ad(Set%File(i)  , TripleCollocation%TC(i))

      call CalcTCforFile_New(Set%File(i) , TripleCollocation%TC(i) , eps_u,eps_v,Scale_u,Scale_v)
      if (Results3C(1:1) /= ' ') write (ResLUN,*) r2_u,eps_u(1),eps_u(2),eps_u(3),r2_v,eps_v(1),eps_v(2),eps_v(3)
      if (Scales3C(1:1) /= ' ') then
        open (ScaleLUN , file=Scales3C , status='unknown')
        write (ScaleLUN,*) Reference,'  ',Scale_u(2),Scale_v(2),Scale_u(3),Scale_v(3)
        close (ScaleLUN , status='keep')
      end if
      if (Errors3C(1:1) /= ' ') then
        open (ErrorLUN , file=Errors3C , status='unknown')
        write (ErrorLUN,*) eps_u(1),eps_v(1),eps_u(2),eps_v(2),eps_u(3),eps_v(3)
        close (ErrorLUN , status='keep')
      end if
      if (Verbosity > 2) call PrintTCMoments
  
      if (.not. StandardOnly) call CalcTCforFileDP_New(Set%File(i) , TripleCollocation%TC(i))
      if (.not. StandardOnly) call CalcTCforFile_Jur(Set%File(i) , TripleCollocation%TC(i))
    enddo
!------------------------------------------------------------------------------!
! File with representation errors as function of s_NWP given:                  !
! apply standard method for all representation errors on first data set.       !
!------------------------------------------------------------------------------!
  else
    open (RepLUN , file=RepresentationErrorFile , status='old')

    read (RepLUN,'(a256)')  In1
    read (RepLUN,'(a256)')  In2
    read (RepLUN,'(a3)')    Par1
    read (RepLUN,'(a3)')    Par2
    read (RepLUN,'(f12.6)') Resol
    read (RepLUN,'(2i8)')   N,Nfour

    write (*,*) ' '
    write (*,*) 'Data from file with representation errors'
    write (*,*) 'File with representation errors  : ',trim(RepresentationErrorFile)
    write (*,*) 'Data file 1 (input for SCORE)    : ',trim(In1)
    write (*,*) 'Data file 2 (input for SCORE)    : ',Trim(In2)
    write (*,*) 'Parameter 1                      : ',Par1
    write (*,*) 'Parameter 2                      : ',Par2
    write (*,*) 'Resolution in km                 :',Resol
    write (*,*) 'Nr of Representation errors      :',N
    write (*,*) 'FFT size of spectral calculation :',Nfour

    do i=1,N
      read (RepLUN,'(5f12.6)') Scale,r2_u,r2_v,r2_t,r2_l
      call CalcTCforFile_New(Set%File(1) , TripleCollocation%TC(1) , eps_u,eps_v,Scale_u,Scale_v)
      if (Results3C(1:1) /= ' ') write (ResLUN,*) Scale, r2_u,eps_u(1),eps_u(2),eps_u(3),r2_v,eps_v(1),eps_v(2),eps_v(3)
    enddo

    close (RepLUN , status='keep')
  end if
!------------------------------------------------------------------------------!
! Finalise.                                                                    !
!------------------------------------------------------------------------------!
  if (Results3C(1:1) /= ' ') close (ResLUN , status='keep')

  end subroutine CalcTripleCollocation


  subroutine CalcTCforFile_Ad(File,TC)
!------------------------------------------------------------------------------!
! Calculates for a single file the moments for triple collocation using only   !
! those points that are not too far from the calibration line.                 !
! When run without bias correction, this version gives the same results as the !
! program of Ad Stoffelen without higher order calibration.                    !
!------------------------------------------------------------------------------!
  implicit none

  type(DataFile_Type), intent(inout)    :: File                     ! Data file struct
  type(Collocation_Type), intent(inout) :: TC                       ! Triple collocation moments struct
!------------------------------------------------------------------------------!
! Local declarations.                                                          !
!------------------------------------------------------------------------------!
  real, dimension(3)                    :: uin,vin                  ! Vector of input (u,v) values
  real, dimension(3)                    :: u,v                      ! Vector of calibrated (u,v) values
  real, dimension(3)                    :: Bias_u,Bias_v            ! Calibration coefficients
  real, dimension(3)                    :: Scale_u,Scale_v          ! Calibration coefficients
  real, dimension(3)                    :: BiasCorrection_u
  real, dimension(3)                    :: BiasCorrection_v
  real, dimension(3)                    :: ScaleCorrection_u
  real, dimension(3)                    :: ScaleCorrection_v
  real, dimension(3)                    :: d2u,d2v                  ! Distance between wind components
  type(Stat_Type),dimension(3)          :: Dist2_u,Dist2_v          ! Distance statistics struct
  real, dimension(3)                    :: d2,d2max                 ! Actual/maximum distance from calibration curve
  real, dimension(3)                    :: eps_u,eps_v              ! Triple collocation error standard deviation
  real, dimension(3)                    :: deps_u,deps_v            ! Precision in triple collocation error standard deviation
  real                                  :: t2_u,t2_v                ! Truth variances
  real                                  :: average,var,abs,min,max  ! Statistics parameters
  real                                  :: RN
  integer                               :: N                        ! Number of points in this file
  integer                               :: Iteration                ! Iteration number
  integer                               :: Rejected                 ! Number of data lines rejected by sigma-test
  integer                               :: j1,j2,j3,l               ! Indices
!------------------------------------------------------------------------------!
! Initialize. Start iteration.                                                 !
!------------------------------------------------------------------------------!
  do j1=1,3
    Bias_u(j1)=0.0
    Bias_v(j1)=0.0

    BiasCorrection_u(j1)=0.0
    BiasCorrection_v(j1)=0.0

    Scale_u(j1)=1.0
    Scale_v(j1)=1.0

    ScaleCorrection_u(j1)=1.0
    ScaleCorrection_v(j1)=1.0

    d2max(j1)=TCfactor**2 * 9.0
  enddo

  write (*,*) ' '
  write (*,*) 'Triple Collocation intermediate results:'
  write (*,*) 'Method          : Ad Stoffelen'
  write (*,*) 'Bias correction :',DoBiasCorrection
  write (*,*) 'Representation error variances in u and v:',r2_u,r2_v

  if (Verbosity > 1) then
    write (*,*) ' '
    select case (Reference)
    case ('buoy')
      write (*,*) '                 B U O Y             S C A T             B A C K'         &
                    //'             S C A T             B A C K'
    case ('scat')
      write (*,*) '                 S C A T             B A C K             B U O Y'         &
                    //'             B A C K             B U O Y'
    case ('back')
      write (*,*) '                 B A C K             B U O Y             S C A T'         &
                    //'             B U O Y             S C A T'
    case default
    end select

    write (*,*) '        Error Std.Dev. 1    Error Std.Dev. 2    Error Std.Dev. 3'         &
                  //'   Scale parameter 2   Scale parameter 3      Truth Std.Dev.'
    write (*,*) 'Iter         u         v         u         v         u         v'         &
                  //'         u         v         u         v         u         v  Rejected'
    write (*,*) '----------------------------------------------------------------'         &
                  //'----------------------------------------------------------------------'
  end if

  Iteration_Loop: do Iteration=1,MaxIter3
!------------------------------------------------------------------------------!
!   Reset first and second moments. Start loop over data points.               !
!------------------------------------------------------------------------------!
    call Reset_Moments(TC%u)
    call Reset_Moments(TC%v)
 
    do j1=1,3
      call Init_Stat(Dist2_u(j1))
      call Init_Stat(Dist2_v(j1))
    enddo

    Rejected=0

    Data_Loop: do l=1,File%NrOfDataLines
      call GetTripleCollocationWinds(File%Line(l) , uin,vin)
!------------------------------------------------------------------------------!
!     Calibration. BiasCorrection remains zero if -nobias option is set.       !
!------------------------------------------------------------------------------!
      do j1=1,3
        u(j1)=(uin(j1) - BiasCorrection_u(j1))*ScaleCorrection_u(j1)
        v(j1)=(vin(j1) - BiasCorrection_v(j1))*ScaleCorrection_v(j1)
      enddo
!------------------------------------------------------------------------------!
!     Sigma test.                                                              !
!------------------------------------------------------------------------------!
      do j1=1,3
        j2=mod(j1,3) + 1
        
        d2u(j1)=(u(j1) - u(j2))**2
        d2v(j1)=(v(j1) - v(j2))**2

        d2(j1)=d2u(j1) + d2v(j1)
        if (d2(j1) >= d2max(j1)) then
          Rejected=Rejected + 1
          cycle Data_Loop
        end if
      enddo
!------------------------------------------------------------------------------!
!     Update first and second moments; update distances                        !
!------------------------------------------------------------------------------!
      call Update_Moments(TC%u , u)
      call Update_Moments(TC%v , v)

      do j1=1,3
        call Update_Stat(Dist2_u(j1),d2u(j1))
        call Update_Stat(Dist2_v(j1),d2v(j1))
      enddo
    enddo Data_Loop
!------------------------------------------------------------------------------!
!   Calculate triple collocation errors.                                       !
!------------------------------------------------------------------------------!
    do j1=1,3
      j2=mod(j1,3) + 1
      j3=mod(j2,3) + 1
      if (j3 == BackIndex) then
        TC%u%M2(1,2)=TC%u%M2(1,2) - r2_u
        TC%v%M2(1,2)=TC%v%M2(1,2) - r2_v
      end if
    enddo

    t2_u=TC%u%M2(1,2)*TC%u%M2(3,1)/TC%u%M2(2,3)
    t2_v=TC%v%M2(1,2)*TC%v%M2(3,1)/TC%v%M2(2,3)

    do j1=1,3
      j2=mod(j1,3) + 1
      eps_u(j1)=sqrt(TC%u%M2(j1,j1) - TC%u%M2(j1,j2))
      eps_v(j1)=sqrt(TC%v%M2(j1,j1) - TC%v%M2(j1,j2))
    enddo

    if (Verbosity > 1) write (*,'(i5,12f10.5,i10)') Iteration,eps_u(1),eps_v(1),eps_u(2),eps_v(2),eps_u(3),eps_v(3), &
                                   Scale_u(2),Scale_v(2),Scale_u(3),Scale_v(3),sqrt(t2_u),sqrt(t2_v),Rejected
!------------------------------------------------------------------------------!
!   Calculate precision in triple collocation error parameters.                !
!------------------------------------------------------------------------------!
    if (Iteration == MaxIter3) then
      RN=1.0/sqrt(2.0*real(TC%u%NrOfPoints))

      do j1=1,3
        j2=mod(j1,3) + 1
        deps_u(j1)=RN*(TC%u%M2(j1,j1) + TC%u%M2(j1,j2))/eps_u(j1)
        deps_v(j1)=RN*(TC%v%M2(j1,j1) + TC%v%M2(j1,j2))/eps_v(j1)
      enddo

     if (Verbosity > 1) write (*,'(a5,6f10.5)') ' Prec',deps_u(1),deps_v(1),deps_u(2),deps_v(2),deps_u(3),deps_v(3)
    end if
!------------------------------------------------------------------------------!
!   Prepare next iteration.                                                    !
!------------------------------------------------------------------------------!
    if (DoBiasCorrection) then
      do j1=1,3
        Bias_u(j1)=TC%u%M1(j1)
        Bias_v(j1)=TC%v%M1(j1)
      enddo
    end if
  
    Scale_u(1)=1.0
    Scale_v(1)=1.0
    Scale_u(2)=TC%u%M2(2,3)/TC%u%M2(3,1)
    Scale_v(2)=TC%v%M2(2,3)/TC%v%M2(3,1)
    Scale_u(3)=TC%u%M2(2,3)/TC%u%M2(1,2)
    Scale_v(3)=TC%v%M2(2,3)/TC%v%M2(1,2)
  
    do j1=1,3
      if (DoBiasCorrection) then
        BiasCorrection_u(j1)=BiasCorrection_u(j1) + Bias_u(j1)
        BiasCorrection_v(j1)=BiasCorrection_v(j1) + Bias_v(j1)
      end if
  
      ScaleCorrection_u(j1)=ScaleCorrection_u(j1)/Scale_u(j1)
      ScaleCorrection_v(j1)=ScaleCorrection_v(j1)/Scale_v(j1)
  
      call Get_Stat(Dist2_u(j1) , average,var,abs,min,max,N , Warning=.false.)
      d2max(j1)=TCfactor**2 * average
      call Get_Stat(Dist2_v(j1) , average,var,abs,min,max,N , Warning=.false.)
      d2max(j1)=d2max(j1) + TCfactor**2 * average
    enddo
  enddo Iteration_Loop
!------------------------------------------------------------------------------!
! Write output.                                                                !
!------------------------------------------------------------------------------!
  write (*,*) ' '
  write (*,*) '           Z O N A L   W I N D   C O M P   U       M E R I D   W I N D   C O M P   V'
  write (*,*) 'Inst      bias     scale     sigma    dsigma      bias     scale     sigma    dsigma'
  write (*,*) '------------------------------------------------------------------------------------'

  select case (Reference)
  case ('buoy')
    write (*,'(a5,8f10.6)') ' Buoy',BiasCorrection_u(1),ScaleCorrection_u(1),eps_u(1),deps_u(1),  &
                                    BiasCorrection_v(1),ScaleCorrection_v(1),eps_v(1),deps_v(1)
    write (*,'(a5,8f10.6)') ' Scat',BiasCorrection_u(2),ScaleCorrection_u(2),eps_u(2),deps_u(2),  &
                                    BiasCorrection_v(2),ScaleCorrection_v(2),eps_v(2),deps_v(2)
    write (*,'(a5,8f10.6)') ' Back',BiasCorrection_u(3),ScaleCorrection_u(3),eps_u(3),deps_u(3),  &
                                    BiasCorrection_v(3),ScaleCorrection_v(3),eps_v(3),deps_v(3)
  case ('scat')
    write (*,'(a5,8f10.6)') ' Buoy',BiasCorrection_u(3),ScaleCorrection_u(3),eps_u(3),deps_u(3),  &
                                    BiasCorrection_v(3),ScaleCorrection_v(3),eps_v(3),deps_v(3)
    write (*,'(a5,8f10.6)') ' Scat',BiasCorrection_u(1),ScaleCorrection_u(1),eps_u(1),deps_u(1),  &
                                    BiasCorrection_v(1),ScaleCorrection_v(1),eps_v(1),deps_v(1)
    write (*,'(a5,8f10.6)') ' Back',BiasCorrection_u(2),ScaleCorrection_u(2),eps_u(2),deps_u(2),  &
                                    BiasCorrection_v(2),ScaleCorrection_v(2),eps_v(2),deps_v(2)
  case ('back')
    write (*,'(a5,8f10.6)') ' Buoy',BiasCorrection_u(2),ScaleCorrection_u(2),eps_u(2),deps_u(2),  &
                                    BiasCorrection_v(2),ScaleCorrection_v(2),eps_v(2),deps_v(2)
    write (*,'(a5,8f10.6)') ' Scat',BiasCorrection_u(3),ScaleCorrection_u(3),eps_u(3),deps_u(3),  &
                                    BiasCorrection_v(3),ScaleCorrection_v(3),eps_v(3),deps_v(3)
    write (*,'(a5,8f10.6)') ' Back',BiasCorrection_u(1),ScaleCorrection_u(1),eps_u(1),deps_u(1),  &
                                    BiasCorrection_v(1),ScaleCorrection_v(1),eps_v(1),deps_v(1)
  case default
  end select

  end subroutine CalcTCforFile_Ad


  subroutine CalcTCforFile_New(File,TC , eps_u,eps_v,ScaleCorrection_u,ScaleCorrection_v)
!------------------------------------------------------------------------------!
! Calculates for a single file the moments for triple collocation using only   !
! those points that are not too far from the calibration line.                 !
! This is  Ad stoffelen's formulation with optimized statistics WHICH IS       !
! CONSIDERED HERE THE STANDARD TRIPLE COLLOCATION METHOD.                      !
!------------------------------------------------------------------------------!
  implicit none

  type(DataFile_Type), intent(inout)    :: File                     ! Data file struct
  type(Collocation_Type), intent(inout) :: TC                       ! Triple collocation moments struct
  real, dimension(3), intent(out)       :: eps_u                    ! Triple collocation error standard deviation for u
  real, dimension(3), intent(out)       :: eps_v                    ! Triple collocation error standard deviation for v
  real, dimension(3), intent(out)       :: ScaleCorrection_u        ! Calibration coefficients (scalings) for u
  real, dimension(3), intent(out)       :: ScaleCorrection_v        ! Calibration coefficients (scalings) for v 
!------------------------------------------------------------------------------!
! Local declarations.                                                          !
!------------------------------------------------------------------------------!
  real, dimension(3)                    :: Dij_u,Dij_v              ! Mixed second order moments
  real, dimension(3)                    :: uin,vin                  ! Vector of input (u,v) values
  real, dimension(3)                    :: u,v                      ! Vector of calibrated (u,v) values
  real, dimension(3)                    :: Bias_u,Bias_v            ! Calibration coefficients (become 0 when iteration converges)
  real, dimension(3)                    :: Scale_u,Scale_v          ! Calibration coefficients (become 1 when iteration converges)
  real, dimension(3)                    :: BiasCorrection_u
  real, dimension(3)                    :: BiasCorrection_v
  real, dimension(3)                    :: dBias_u,dBias_v
  real, dimension(3)                    :: dScale_u,dScale_v
  real, dimension(3)                    :: d2u,d2v                  ! Distance between wind components
  type(Stat_Type),dimension(3)          :: Dist2_u,Dist2_v          ! Distance statistics struct
  real, dimension(3)                    :: d2,d2max                 ! Actual/maximum distance from calibration curve
  real, dimension(3)                    :: deps_u,deps_v            ! Precision in triple collocation error standard deviation
  real                                  :: t2_u,t2_v                ! Truth variances
  real                                  :: average,var,abe,min,max  ! Statistics parameters
  real                                  :: RN
  real                                  :: AverageAngle,Angle
  integer                               :: N                        ! Number of points in this file
  integer                               :: Iteration                ! Iteration number
  integer                               :: Rejected                 ! Number of data lines rejected by sigma-test
  integer                               :: j1,j2,j3,l               ! Indices
!------------------------------------------------------------------------------!
! Initialize. Start iteration.                                                 !
!------------------------------------------------------------------------------!
  do j1=1,3
    Bias_u(j1)=0.0
    Bias_v(j1)=0.0

    BiasCorrection_u(j1)=0.0
    BiasCorrection_v(j1)=0.0

    Scale_u(j1)=1.0
    Scale_v(j1)=1.0

    ScaleCorrection_u(j1)=1.0
    ScaleCorrection_v(j1)=1.0

    d2max(j1)=TCfactor**2 * 9.0
  enddo

  write (*,*) ' '
  write (*,*) 'Triple Collocation results:'
  write (*,*) 'Method          : Ad Stoffelen with optimized statistics (New)'
  write (*,*) 'Bias correction :',DoBiasCorrection
  write (*,'(a,2f9.6)') ' Representation error variances      in u and v :',r2_u,r2_v
  write (*,'(a,2f9.6)') ' Scatterometer white noise variances in u and v :',wnv_u,wnv_v
  write (*,'(a,2f9.6)') ' Correlations between scat and back  in u and v :',csb_u,csb_v
  write (*,*) 'Use only collocations common with other data set?',Common3C

  if (Verbosity > 1) then
    write (*,*) ' '
    select case (Reference)
    case ('buoy')
      write (*,*) '                 B U O Y             S C A T             B A C K'         &
                    //'             S C A T             B A C K'
    case ('scat')
      write (*,*) '                 S C A T             B A C K             B U O Y'         &
                    //'             B A C K             B U O Y'
    case ('back')
      write (*,*) '                 B A C K             B U O Y             S C A T'         &
                    //'             B U O Y             S C A T'
    case default
    end select

    write (*,*) '        Error Std.Dev. 1    Error Std.Dev. 2    Error Std.Dev. 3'         &
                  //'   Scale parameter 2   Scale parameter 3      Truth Std.Dev.'
    write (*,*) 'Iter         u         v         u         v         u         v'         &
                  //'         u         v         u         v         u         v  Rejected'
    write (*,*) '----------------------------------------------------------------'         &
                  //'----------------------------------------------------------------------'
  end if

  Iteration_Loop: do Iteration=1,MaxIter3
!------------------------------------------------------------------------------!
!   Reset first and second moments. Start loop over data points.               !
!------------------------------------------------------------------------------!
    call Reset_Moments(TC%u)
    call Reset_Moments(TC%v)

    do j1=1,3
      Dij_u(j1)=0.0
      Dij_v(j1)=0.0
      call Init_Stat(Dist2_u(j1))
      call Init_Stat(Dist2_v(j1))
    enddo

    Rejected=0
    N=0

    Data_Loop: do l=1,File%NrOfDataLines
      if (Common3C .and. .not. File%Line(l)%Common) cycle                      ! Reject collocations that are not common if Common3C
      call GetTripleCollocationWinds(File%Line(l) , uin,vin)

      AverageAngle=0.5*(File%Line(l)%BuoyDir + File%Line(l)%BackDir)
      Angle=abs(AverageAngle - File%Line(l)%ScatDir)
      if (MaxAngle > 0 .and. Angle > MaxAngle) cycle
!------------------------------------------------------------------------------!
!     Calibration. BiasCorrection remains zero if -nobias option is set.       !
!------------------------------------------------------------------------------!
      do j1=1,3
        u(j1)=(uin(j1) - BiasCorrection_u(j1))*ScaleCorrection_u(j1)
        v(j1)=(vin(j1) - BiasCorrection_v(j1))*ScaleCorrection_v(j1)
      enddo
!------------------------------------------------------------------------------!
!     Sigma test.                                                              !
!------------------------------------------------------------------------------!
      do j1=1,3
        j2=mod(j1,3) + 1

        d2u(j1)=(u(j1) - u(j2))**2
        d2v(j1)=(v(j1) - v(j2))**2

        d2(j1)=d2u(j1) + d2v(j1)
        if (d2(j1) >= d2max(j1)) then
          Rejected=Rejected + 1
          cycle Data_Loop
        end if
      enddo
!------------------------------------------------------------------------------!
!     Update optimized second moments and distances                            !
!------------------------------------------------------------------------------!
      call Update_Moments(TC%u , u)
      call Update_Moments(TC%v , v)

      N=N + 1
      RN=1.0/real(N)

      do j1=1,3
        j2=mod(j1,3) + 1
       
        Dij_u(j1)=Dij_u(j1) + RN*(u(j1)*(u(j1) - u(j2)) - Dij_u(j1))
        Dij_v(j1)=Dij_v(j1) + RN*(v(j1)*(v(j1) - v(j2)) - Dij_v(j1))
        call Update_Stat(Dist2_u(j1),d2u(j1))
        call Update_Stat(Dist2_v(j1),d2v(j1))
      enddo
    enddo Data_Loop
!------------------------------------------------------------------------------!
!   Calculate triple collocation errors.                                       !
!------------------------------------------------------------------------------!
    do j1=1,3
      j2=mod(j1,3) + 1
      j3=mod(j2,3) + 1

      if (j3 == BackIndex) then                   ! then j1 == BuoyIndex and j2 == ScatIndex;
        Dij_u(j1)=Dij_u(j1) - r2_u                ! correct for representation error
        Dij_v(j1)=Dij_v(j1) - r2_v
        TC%u%M2(j1,j2)=TC%u%M2(j1,j2) - r2_u
        TC%v%M2(j1,j2)=TC%v%M2(j1,j2) - r2_v
      end if

      if (j3 == BuoyIndex) then                   ! then j1 == ScatIndex and j2 == BackIndex;
        Dij_u(j1)=Dij_u(j1) - csb_u               ! correct for correlation between scat and buoy
        Dij_v(j1)=Dij_v(j1) - csb_v
        TC%u%M2(j1,j2)=TC%u%M2(j1,j2) - csb_u
        TC%v%M2(j1,j2)=TC%v%M2(j1,j2) - csb_v
      end if

      if (j1 == ScatIndex) then
        Dij_u(j1)=Dij_u(j1) - wnv_u               ! Correct for white noise contribution in scat
        Dij_v(j1)=Dij_v(j1) - wnv_v
        TC%u%M2(j1,j1)=TC%u%M2(j1,j1) - wnv_u
        TC%v%M2(j1,j1)=TC%v%M2(j1,j1) - wnv_v
      end if
    enddo

    t2_u=TC%u%M2(1,2)*TC%u%M2(3,1)/TC%u%M2(2,3)
    t2_v=TC%v%M2(1,2)*TC%v%M2(3,1)/TC%v%M2(2,3)

    do j1=1,3
      j2=mod(j1,3) + 1
!!    eps_u(j1)=sqrt(Dij_u(j1))                   ! Om een of andere reden gaat dit fout
!!    eps_v(j1)=sqrt(Dij_v(j1))                   ! voor j1 = 1 bij referentie buoy

      eps_u(j1)=sqrt(TC%u%M2(j1,j1) - TC%u%M2(j1,j2))
      eps_v(j1)=sqrt(TC%v%M2(j1,j1) - TC%v%M2(j1,j2))
    enddo

    if (Verbosity > 1) write (*,'(i5,12f10.5,i10)') Iteration,eps_u(1),eps_v(1),eps_u(2),eps_v(2),eps_u(3),eps_v(3), &
                                   Scale_u(2),Scale_v(2),Scale_u(3),Scale_v(3),sqrt(t2_u),sqrt(t2_v),Rejected
!------------------------------------------------------------------------------!
!   Calculate precision in triple collocation error parameters.                !
!------------------------------------------------------------------------------!
    if (Iteration == MaxIter3) then
      RN=1.0/sqrt(real(TC%u%NrOfPoints))

      do j1=1,3
        j2=mod(j1,3) + 1
        deps_u(j1)=RN*eps_u(j1)*sqrt(2.0*eps_u(j1)**2 + eps_u(j2)**2)
        deps_v(j1)=RN*eps_v(j1)*sqrt(2.0*eps_v(j1)**2 + eps_v(j2)**2)
      enddo

      dBias_u(1)=RN*eps_u(1)
      dBias_v(1)=RN*eps_v(1)
      dBias_u(2)=RN*eps_u(2)
      dBias_v(2)=RN*eps_v(2)
      dBias_u(3)=RN*eps_u(3)
      dBias_v(3)=RN*eps_v(3)

      dScale_u(1)=0.0
      dScale_v(1)=0.0
      dScale_u(2)=RN/TC%u%M2(2,3) * sqrt(eps_u(2)**2*eps_u(3)**2 + eps_u(3)**2*eps_u(1)**2 + &  ! Should all be minus, but
                                         6.0*eps_u(1)**2*eps_u(2)**2*eps_u(3)**4)               ! that results in NaN's
      dScale_v(2)=RN/TC%v%M2(2,3) * sqrt(eps_v(2)**2*eps_v(3)**2 + eps_v(3)**2*eps_v(1)**2 + &
                                         6.0*eps_v(1)**2*eps_v(2)**2*eps_v(3)**4)
      dScale_u(3)=RN/TC%u%M2(2,3) * sqrt(eps_u(2)**2*eps_u(3)**2 + eps_u(1)**2*eps_u(2)**2 + &
                                         6.0*eps_u(1)**2*eps_u(2)**4*eps_u(3)**2)
      dScale_v(3)=RN/TC%v%M2(2,3) * sqrt(eps_v(2)**2*eps_v(3)**2 + eps_v(1)**2*eps_v(2)**2 + &
                                         6.0*eps_v(1)**2*eps_v(2)**4*eps_v(3)**2)

      if (Verbosity > 1) write (*,'(a5,10f10.5)') ' Prec',deps_u(1),deps_v(1),deps_u(2),deps_v(2),deps_u(3),deps_v(3), &
                                                    dScale_u(2),dScale_v(2),dScale_u(3),dScale_v(3)
    end if
!------------------------------------------------------------------------------!
!   Prepare next iteration.                                                    !
!------------------------------------------------------------------------------!
    if (DoBiasCorrection) then
      do j1=1,3
        Bias_u(j1)=TC%u%M1(j1)
        Bias_v(j1)=TC%v%M1(j1)
      enddo
    end if

    Scale_u(1)=1.0
    Scale_v(1)=1.0
    Scale_u(2)=TC%u%M2(2,3)/TC%u%M2(3,1)
    Scale_v(2)=TC%v%M2(2,3)/TC%v%M2(3,1)
    Scale_u(3)=TC%u%M2(2,3)/TC%u%M2(1,2)
    Scale_v(3)=TC%v%M2(2,3)/TC%v%M2(1,2)

    do j1=1,3
      if (DoBiasCorrection) then
        BiasCorrection_u(j1)=BiasCorrection_u(j1) + Bias_u(j1)
        BiasCorrection_v(j1)=BiasCorrection_v(j1) + Bias_v(j1)
      end if

      ScaleCorrection_u(j1)=ScaleCorrection_u(j1)/Scale_u(j1)
      ScaleCorrection_v(j1)=ScaleCorrection_v(j1)/Scale_v(j1)

      call Get_Stat(Dist2_u(j1) , average,var,abe,min,max,N , Warning=.false.)
      d2max(j1)=TCfactor**2 * average
      call Get_Stat(Dist2_v(j1) , average,var,abe,min,max,N , Warning=.false.)
      d2max(j1)=d2max(j1) + TCfactor**2 * average
    enddo
  enddo Iteration_Loop
!------------------------------------------------------------------------------!
! Write output.                                                                !
!------------------------------------------------------------------------------!
  write (*,*) ' '
  write (*,*) '                               Z O N A L   W I N D   C O M P   U' //  &
                  '                           M E R I D   W I N D   C O M P   V'
  write (*,*) 'Inst      bias     dbias     scale    dscale     sigma    dsigma' //  &
                  '      bias     dbias     scale    dscale     sigma    dsigma'
  write (*,*) '----------------------------------------------------------------' //   &
                  '------------------------------------------------------------'

  select case (Reference)
  case ('buoy')
    write (*,'(a5,12f10.6)') ' Buoy',BiasCorrection_u(1),dBias_u(1),ScaleCorrection_u(1),dScale_u(1),eps_u(1),deps_u(1),  &
                                     BiasCorrection_v(1),dBias_v(1),ScaleCorrection_v(1),dScale_v(1),eps_v(1),deps_v(1)
    write (*,'(a5,12f10.6)') ' Scat',BiasCorrection_u(2),dBias_u(2),ScaleCorrection_u(2),dScale_u(2),eps_u(2),deps_u(2),  &
                                     BiasCorrection_v(2),dBias_v(2),ScaleCorrection_v(2),dScale_v(2),eps_v(2),deps_v(2)
    write (*,'(a5,12f10.6)') ' Back',BiasCorrection_u(3),dBias_u(3),ScaleCorrection_u(3),dScale_u(3),eps_u(3),deps_u(3),  &
                                     BiasCorrection_v(3),dBias_v(3),ScaleCorrection_v(3),dScale_v(3),eps_v(3),deps_v(3)
  case ('scat')
    write (*,'(a5,12f10.6)') ' Buoy',BiasCorrection_u(3),dBias_u(3),ScaleCorrection_u(3),dScale_u(3),eps_u(3),deps_u(3),  &
                                     BiasCorrection_v(3),dBias_v(3),ScaleCorrection_v(3),dScale_v(3),eps_v(3),deps_v(3)
    write (*,'(a5,12f10.6)') ' Scat',BiasCorrection_u(1),dBias_u(1),ScaleCorrection_u(1),dScale_u(1),eps_u(1),deps_u(1),  &
                                     BiasCorrection_v(1),dBias_v(1),ScaleCorrection_v(1),dScale_v(1),eps_v(1),deps_v(1)
    write (*,'(a5,12f10.6)') ' Back',BiasCorrection_u(2),dBias_u(2),ScaleCorrection_u(2),dScale_u(2),eps_u(2),deps_u(2),  &
                                     BiasCorrection_v(2),dBias_v(2),ScaleCorrection_v(2),dScale_v(2),eps_v(2),deps_v(2)
  case ('back')
    write (*,'(a5,12f10.6)') ' Buoy',BiasCorrection_u(2),dBias_u(2),ScaleCorrection_u(2),dScale_u(2),eps_u(2),deps_u(2),  &
                                     BiasCorrection_v(2),dBias_v(2),ScaleCorrection_v(2),dScale_v(2),eps_v(2),deps_v(2)
    write (*,'(a5,12f10.6)') ' Scat',BiasCorrection_u(3),dBias_u(3),ScaleCorrection_u(3),dScale_u(3),eps_u(3),deps_u(3),  &
                                     BiasCorrection_v(3),dBias_v(3),ScaleCorrection_v(3),dScale_v(3),eps_v(3),deps_v(3)
    write (*,'(a5,12f10.6)') ' Back',BiasCorrection_u(1),dBias_u(1),ScaleCorrection_u(1),dScale_u(1),eps_u(1),deps_u(1),  &
                                     BiasCorrection_v(1),dBias_v(1),ScaleCorrection_v(1),dScale_v(1),eps_v(1),deps_v(1)
  case default
  end select

  end subroutine CalcTCforFile_New


  subroutine CalcTCforFileDP_New(File,TC)
!------------------------------------------------------------------------------!
! Calculates for a single file the moments for triple collocation using only   !
! those points that are not too far from the calibration line.                 !
! This is  Ad stoffelen's formulation with optimized statistics in double      !
! precision.                                                                   !
!------------------------------------------------------------------------------!
  implicit none

  type(DataFile_Type), intent(inout)    :: File                     ! Data file struct
  type(Collocation_Type), intent(inout) :: TC                       ! Triple collocation moments struct
!------------------------------------------------------------------------------!
! Local declarations.                                                          !
!------------------------------------------------------------------------------!
  double precision, dimension(3)        :: Dij_u,Dij_v              ! Mixed second order moments
  double precision                      :: DN
  real, dimension(3)                    :: uin,vin                  ! Vector of input (u,v) values
  real, dimension(3)                    :: u,v                      ! Vector of calibrated (u,v) values
  real, dimension(3)                    :: Bias_u,Bias_v            ! Calibration coefficients
  real, dimension(3)                    :: Scale_u,Scale_v          ! Calibration coefficients
  real, dimension(3)                    :: BiasCorrection_u
  real, dimension(3)                    :: BiasCorrection_v
  real, dimension(3)                    :: ScaleCorrection_u
  real, dimension(3)                    :: ScaleCorrection_v
  real, dimension(3)                    :: d2u,d2v                  ! Distance between wind components
  type(Stat_Type),dimension(3)          :: Dist2_u,Dist2_v          ! Distance statistics struct
  real, dimension(3)                    :: d2,d2max                 ! Actual/maximum distance from calibration curve
  real, dimension(3)                    :: eps_u,eps_v              ! Triple collocation error standard deviation
  real, dimension(3)                    :: deps_u,deps_v            ! Precision in triple collocation error standard deviation
  real                                  :: t2_u,t2_v                ! Truth variances
  real                                  :: average,var,abs,min,max  ! Statistics parameters
  real                                  :: RN
  integer                               :: N                        ! Number of points in this file
  integer                               :: Iteration                ! Iteration number
  integer                               :: Rejected                 ! Number of data lines rejected by sigma-test
  integer                               :: j1,j2,j3,l               ! Indices
!------------------------------------------------------------------------------!
! Initialize. Start iteration.                                                 !
!------------------------------------------------------------------------------!
  do j1=1,3
    Bias_u(j1)=0.0
    Bias_v(j1)=0.0

    BiasCorrection_u(j1)=0.0
    BiasCorrection_v(j1)=0.0

    Scale_u(j1)=1.0
    Scale_v(j1)=1.0

    ScaleCorrection_u(j1)=1.0
    ScaleCorrection_v(j1)=1.0

    d2max(j1)=TCfactor**2 * 9.0
  enddo

  write (*,*) ' '
  write (*,*) 'Triple Collocation results:'
  write (*,*) 'Method          : Ad Stoffelen with optimized statistics in double precision'
  write (*,*) 'Bias correction :',DoBiasCorrection
  write (*,*) 'Representation error variances in u and v:',r2_u,r2_v

  if (Verbosity > 1) then
    write (*,*) ' '
    select case (Reference)
    case ('buoy')
      write (*,*) '                 B U O Y             S C A T             B A C K'         &
                    //'             S C A T             B A C K'
    case ('scat')
      write (*,*) '                 S C A T             B A C K             B U O Y'         &
                    //'             B A C K             B U O Y'
    case ('back')
      write (*,*) '                 B A C K             B U O Y             S C A T'         &
                    //'             B U O Y             S C A T'
    case default
    end select

    write (*,*) '        Error Std.Dev. 1    Error Std.Dev. 2    Error Std.Dev. 3'         &
                  //'   Scale parameter 2   Scale parameter 3      Truth Std.Dev.'
    write (*,*) 'Iter         u         v         u         v         u         v'         &
                  //'         u         v         u         v         u         v  Rejected'
    write (*,*) '----------------------------------------------------------------'         &
                  //'----------------------------------------------------------------------'
  end if

  Iteration_Loop: do Iteration=1,MaxIter3
!------------------------------------------------------------------------------!
!   Reset first and second moments. Start loop over data points.               !
!------------------------------------------------------------------------------!
    call Reset_Moments(TC%u)
    call Reset_Moments(TC%v)

    do j1=1,3
      Dij_u(j1)=0.0
      Dij_v(j1)=0.0
      call Init_Stat(Dist2_u(j1))
      call Init_Stat(Dist2_v(j1))
    enddo

    Rejected=0
    N=0

    Data_Loop: do l=1,File%NrOfDataLines
      call GetTripleCollocationWinds(File%Line(l) , uin,vin)
!------------------------------------------------------------------------------!
!     Calibration. BiasCorrection remains zero if -nobias option is set.       !
!------------------------------------------------------------------------------!
      do j1=1,3
        u(j1)=(uin(j1) - BiasCorrection_u(j1))*ScaleCorrection_u(j1)
        v(j1)=(vin(j1) - BiasCorrection_v(j1))*ScaleCorrection_v(j1)
      enddo
!------------------------------------------------------------------------------!
!     Sigma test.                                                              !
!------------------------------------------------------------------------------!
      do j1=1,3
        j2=mod(j1,3) + 1

        d2u(j1)=(u(j1) - u(j2))**2
        d2v(j1)=(v(j1) - v(j2))**2

        d2(j1)=d2u(j1) + d2v(j1)
        if (d2(j1) >= d2max(j1)) then
          Rejected=Rejected + 1
          cycle Data_Loop
        end if
      enddo
!------------------------------------------------------------------------------!
!     Update optimized second moments and distances                            !
!------------------------------------------------------------------------------!
      call Update_Moments(TC%u , u)
      call Update_Moments(TC%v , v)

      N=N + 1
      DN=1.0/dble(N)

      do j1=1,3
        j2=mod(j1,3) + 1

        Dij_u(j1)=Dij_u(j1) + DN*(dble(u(j1)*(u(j1) - u(j2))) - Dij_u(j1))
        Dij_v(j1)=Dij_v(j1) + DN*(dble(v(j1)*(v(j1) - v(j2))) - Dij_v(j1))
        call Update_Stat(Dist2_u(j1),d2u(j1))
        call Update_Stat(Dist2_v(j1),d2v(j1))
      enddo
    enddo Data_Loop
!------------------------------------------------------------------------------!
!   Calculate triple collocation errors.                                       !
!------------------------------------------------------------------------------!
    do j1=1,3
      j2=mod(j1,3) + 1
      j3=mod(j2,3) + 1
      if (j3 == BackIndex) then
        Dij_u(j1)=Dij_u(j1) - dble(r2_u)
        Dij_v(j1)=Dij_v(j1) - dble(r2_v)
        TC%u%M2(j1,j2)=TC%u%M2(j1,j2) - r2_u
        TC%v%M2(j1,j2)=TC%v%M2(j1,j2) - r2_v
      end if
    enddo

    t2_u=TC%u%M2(1,2)*TC%u%M2(3,1)/TC%u%M2(2,3)
    t2_v=TC%v%M2(1,2)*TC%v%M2(3,1)/TC%v%M2(2,3)

    do j1=1,3
      j2=mod(j1,3) + 1
!!    eps_u(j1)=sqrt(real(Dij_u(j1)))                  ! Gaat fout voor j1=1 bij referentie buoy
!!    eps_v(j1)=sqrt(real(Dij_v(j1)))

      eps_u(j1)=sqrt(TC%u%M2(j1,j1) - TC%u%M2(j1,j2))
      eps_v(j1)=sqrt(TC%v%M2(j1,j1) - TC%v%M2(j1,j2))
    enddo

    if (Verbosity > 1) write (*,'(i5,12f10.5,i10)') Iteration,eps_u(1),eps_v(1),eps_u(2),eps_v(2),eps_u(3),eps_v(3), &
                                   Scale_u(2),Scale_v(2),Scale_u(3),Scale_v(3),sqrt(t2_u),sqrt(t2_v),Rejected
!------------------------------------------------------------------------------!
!   Calculate precision in triple collocation error parameters.                !
!------------------------------------------------------------------------------!
    if (Iteration == MaxIter3) then
      RN=1.0/sqrt(real(TC%u%NrOfPoints))

      do j1=1,3
        j2=mod(j1,3) + 1
        deps_u(j1)=RN*eps_u(j1)*sqrt(2.0*eps_u(j1)**2 + eps_u(j2)**2)
        deps_v(j1)=RN*eps_v(j1)*sqrt(2.0*eps_v(j1)**2 + eps_v(j2)**2)
      enddo

      if (Verbosity > 1) write (*,'(a5,6f10.5)') ' Prec',deps_u(1),deps_v(1),deps_u(2),deps_v(2),deps_u(3),deps_v(3)
    end if
!------------------------------------------------------------------------------!
!   Prepare next iteration.                                                    !
!------------------------------------------------------------------------------!
    if (DoBiasCorrection) then
      do j1=1,3
        Bias_u(j1)=TC%u%M1(j1)
        Bias_v(j1)=TC%v%M1(j1)
      enddo
    end if

    Scale_u(1)=1.0
    Scale_v(1)=1.0
    Scale_u(2)=TC%u%M2(2,3)/TC%u%M2(3,1)
    Scale_v(2)=TC%v%M2(2,3)/TC%v%M2(3,1)
    Scale_u(3)=TC%u%M2(2,3)/TC%u%M2(1,2)
    Scale_v(3)=TC%v%M2(2,3)/TC%v%M2(1,2)

    do j1=1,3
      if (DoBiasCorrection) then
        BiasCorrection_u(j1)=BiasCorrection_u(j1) + Bias_u(j1)
        BiasCorrection_v(j1)=BiasCorrection_v(j1) + Bias_v(j1)
      end if

      ScaleCorrection_u(j1)=ScaleCorrection_u(j1)/Scale_u(j1)
      ScaleCorrection_v(j1)=ScaleCorrection_v(j1)/Scale_v(j1)

      call Get_Stat(Dist2_u(j1) , average,var,abs,min,max,N , Warning=.false.)
      d2max(j1)=TCfactor**2 * average
      call Get_Stat(Dist2_v(j1) , average,var,abs,min,max,N , Warning=.false.)
      d2max(j1)=d2max(j1) + TCfactor**2 * average
    enddo
  enddo Iteration_Loop
!------------------------------------------------------------------------------!
! Write output.                                                                !
!------------------------------------------------------------------------------!
  write (*,*) ' '
  write (*,*) '           Z O N A L   W I N D   C O M P   U       M E R I D   W I N D   C O M P   V'
  write (*,*) 'Inst      bias     scale     sigma    dsigma      bias     scale     sigma    dsigma'
  write (*,*) '------------------------------------------------------------------------------------'

  select case (Reference)
  case ('buoy')
    write (*,'(a5,8f10.6)') ' Buoy',BiasCorrection_u(1),ScaleCorrection_u(1),eps_u(1),deps_u(1),  &
                                    BiasCorrection_v(1),ScaleCorrection_v(1),eps_v(1),deps_v(1)
    write (*,'(a5,8f10.6)') ' Scat',BiasCorrection_u(2),ScaleCorrection_u(2),eps_u(2),deps_u(2),  &
                                    BiasCorrection_v(2),ScaleCorrection_v(2),eps_v(2),deps_v(2)
    write (*,'(a5,8f10.6)') ' Back',BiasCorrection_u(3),ScaleCorrection_u(3),eps_u(3),deps_u(3),  &
                                    BiasCorrection_v(3),ScaleCorrection_v(3),eps_v(3),deps_v(3)
  case ('scat')
    write (*,'(a5,8f10.6)') ' Buoy',BiasCorrection_u(3),ScaleCorrection_u(3),eps_u(3),deps_u(3),  &
                                    BiasCorrection_v(3),ScaleCorrection_v(3),eps_v(3),deps_v(3)
    write (*,'(a5,8f10.6)') ' Scat',BiasCorrection_u(1),ScaleCorrection_u(1),eps_u(1),deps_u(1),  &
                                    BiasCorrection_v(1),ScaleCorrection_v(1),eps_v(1),deps_v(1)
    write (*,'(a5,8f10.6)') ' Back',BiasCorrection_u(2),ScaleCorrection_u(2),eps_u(2),deps_u(2),  &
                                    BiasCorrection_v(2),ScaleCorrection_v(2),eps_v(2),deps_v(2)
  case ('back')
    write (*,'(a5,8f10.6)') ' Buoy',BiasCorrection_u(2),ScaleCorrection_u(2),eps_u(2),deps_u(2),  &
                                    BiasCorrection_v(2),ScaleCorrection_v(2),eps_v(2),deps_v(2)
    write (*,'(a5,8f10.6)') ' Scat',BiasCorrection_u(3),ScaleCorrection_u(3),eps_u(3),deps_u(3),  &
                                    BiasCorrection_v(3),ScaleCorrection_v(3),eps_v(3),deps_v(3)
    write (*,'(a5,8f10.6)') ' Back',BiasCorrection_u(1),ScaleCorrection_u(1),eps_u(1),deps_u(1),  &
                                    BiasCorrection_v(1),ScaleCorrection_v(1),eps_v(1),deps_v(1)
  case default
  end select

  end subroutine CalcTCforFileDP_New


  subroutine CalcTCforFile_Jur(File,TC)
!------------------------------------------------------------------------------!
! Calculates for a single file the moments for triple collocation using only   !
! those points that are not too far from the calibration line.                 !
! This version uses the formulas of Jur Vogelzang.                             !
!------------------------------------------------------------------------------!
  implicit none

  type(DataFile_Type), intent(inout)    :: File                     ! Data file struct
  type(Collocation_Type), intent(inout) :: TC                       ! Triple collocation moments struct
!------------------------------------------------------------------------------!
! Local declarations.                                                          !
!------------------------------------------------------------------------------!
  real, dimension(3)                    :: u,v                      ! Vector of (u,v) values
  real, dimension(3)                    :: a_u,b_u,a_v,b_v          ! Calibration coefficients
  real, dimension(3)                    :: d2u,d2v                  ! Distance between wind components
  type(Stat_Type),dimension(3)          :: Dist2_u,Dist2_v          ! Distance statistics struct
  real, dimension(3)                    :: d2,d2max                 ! Actual/maximum distance from calibration curve
  real, dimension(3)                    :: eps_u,eps_v              ! Triple collocation errors
  real                                  :: t2_u,t2_v                ! Truth variances
  real, dimension(3)                    :: deps_u,deps_v            ! Precision in triple collocation error estimates
  real                                  :: RN
  real, dimension(3,3)                  :: C_u,C_v                  ! Covariances
  real                                  :: average,var,abs,min,max  ! Statistics parameters
  integer                               :: N                        ! Number of points in this file
  integer                               :: Iteration                ! Iteration number
  integer                               :: Rejected                 ! Number of data lines rejected by sigma-test
  integer                               :: j1,j2,j3,l               ! Indices
!------------------------------------------------------------------------------!
! Initialize. Start iteration.                                                 !
!------------------------------------------------------------------------------!
  do j1=1,3
    b_u(j1)=0.0
    b_v(j1)=0.0

    a_u(j1)=1.0
    a_v(j1)=1.0

    d2max(j1)=TCfactor**2 * 9.0
  enddo

  write (*,*) ' '
  write (*,*) 'Triple Collocation results:'
  write (*,*) 'Method          : Jur Vogelzang'
  write (*,*) 'Bias correction :',DoBiasCorrection
  write (*,*) 'Representation error variances in u and v:',r2_u,r2_v

  if (Verbosity > 1) then
    write (*,*) ' '
    select case (Reference)
    case ('buoy')
      write (*,*) '                 B U O Y             S C A T             B A C K'         &
                    //'             S C A T             B A C K'
    case ('scat')
      write (*,*) '                 S C A T             B A C K             B U O Y'         &
                    //'             B A C K             B U O Y'
    case ('back')
      write (*,*) '                 B A C K             B U O Y             S C A T'         &
                    //'             B U O Y             S C A T'
    case default
    end select

    write (*,*) '        Error Std.Dev. 1    Error Std.Dev. 2    Error Std.Dev. 3'         &
                  //'   Scale parameter 2   Scale parameter 3      Truth Std.Dev.'
    write (*,*) 'Iter         u         v         u         v         u         v'         &
                  //'         u         v         u         v         u         v  Rejected'
    write (*,*) '----------------------------------------------------------------'         &
                  //'----------------------------------------------------------------------'
  end if

  Iteration_Loop: do Iteration=1,MaxIter3
!------------------------------------------------------------------------------!
!   Reset first and second moments. Start loop over data points.               !
!------------------------------------------------------------------------------!
    call Reset_Moments(TC%u)
    call Reset_Moments(TC%v)

    do j1=1,3
      call Init_Stat(Dist2_u(j1))
      call Init_Stat(Dist2_v(j1))
    enddo

    Rejected=0

    Data_Loop: do l=1,File%NrOfDataLines
      call GetTripleCollocationWinds(File%Line(l) , u,v)
!------------------------------------------------------------------------------!
!     Sigma test.                                                              !
!------------------------------------------------------------------------------!
      do j1=1,3
        j2=mod(j1,3) + 1

        d2u(j1)=(b_u(j1) + a_u(j1)*u(j1) - b_u(j2) - a_u(j2)*u(j2))**2
        d2v(j1)=(b_v(j1) + a_v(j1)*v(j1) - b_v(j2) - a_v(j2)*v(j2))**2

        d2(j1)=d2u(j1) + d2v(j1)
        if (d2(j1) >= d2max(j1)) then
          Rejected=Rejected + 1
          cycle Data_Loop
        end if
      enddo
!------------------------------------------------------------------------------!
!     Update first and second moments; update distances                        !
!------------------------------------------------------------------------------!
      call Update_Moments(TC%u , u)
      call Update_Moments(TC%v , v)

      do j1=1,3
        call Update_Stat(Dist2_u(j1),d2u(j1))
        call Update_Stat(Dist2_v(j1),d2v(j1))
      enddo
    enddo Data_Loop
!------------------------------------------------------------------------------!
!   Calculate triple collocation errors.                                       !
!------------------------------------------------------------------------------!
    if (DoBiasCorrection) then
      do j1=1,3
        do j2=1,3
          C_u(j1,j2)=TC%u%M2(j1,j2) - TC%u%M1(j1)*TC%u%M1(j2)  ! Set variances and covariances
          C_v(j1,j2)=TC%v%M2(j1,j2) - TC%v%M1(j1)*TC%v%M1(j2)
        enddo
        j2=mod(j1,3) + 1
        j3=mod(j2,3) + 1
        if (j3 == BackIndex) then
          C_u(j1,j2)=C_u(j1,j2) - r2_u                         ! Correct for representation error
          C_v(j1,j2)=C_v(j1,j2) - r2_v
        end if
      enddo
    else
      do j1=1,3
        do j2=1,3
          C_u(j1,j2)=TC%u%M2(j1,j2)
          C_v(j1,j2)=TC%v%M2(j1,j2)
        enddo
        j2=mod(j1,3) + 1
        j3=mod(j2,3) + 1
        if (j3 == BackIndex) then
          C_u(j1,j2)=C_u(j1,j2) - r2_u
          C_v(j1,j2)=C_v(j1,j2) - r2_v
        end if
      enddo
    end if

    t2_u=C_u(1,2)*C_u(3,1)/C_u(2,3)
    t2_v=C_v(1,2)*C_v(3,1)/C_v(2,3)

    do j1=1,3
      j2=mod(j1,3) + 1
      j3=mod(j2,3) + 1

      b_u(j1)=TC%u%M1(j1)
      b_v(j1)=TC%v%M1(j1)

      a_u(j1)=C_u(2,3)/C_u(j2,j3)
      a_v(j1)=C_v(2,3)/C_v(j2,j3)

      eps_u(j1)=sqrt(C_u(j1,j1) - a_u(j1)**2 * t2_u)
      eps_v(j1)=sqrt(C_v(j1,j1) - a_v(j1)**2 * t2_v)
    enddo

    if (Verbosity > 1) write (*,'(i5,12f10.5,i10)') Iteration,eps_u(1),eps_v(1),eps_u(2),eps_v(2),eps_u(3),eps_v(3), &
                                   a_u(2),a_v(2),a_u(3),a_v(3),sqrt(t2_u),sqrt(t2_v),Rejected
!------------------------------------------------------------------------------!
!   Calculate precision in triple collocation error estimates.                 !
!------------------------------------------------------------------------------!
    if (Iteration == MaxIter3) then
      RN=sqrt(1.0/(2.0*real(TC%u%NrOfPoints)))

      do j1=1,3
        deps_u(j1)=RN*(C_u(j1,j1) + 3.0*t2_u)/eps_u(j1)
        deps_v(j1)=RN*(C_v(j1,j1) + 3.0*t2_v)/eps_v(j1)
      enddo

      if (Verbosity > 1) write (*,'(a5,6f10.5)') ' Prec',deps_u(1),deps_v(1),deps_u(2),deps_v(2),deps_u(3),deps_v(3)
    end if
!------------------------------------------------------------------------------!
!   Prepare next iteration.                                                    !
!------------------------------------------------------------------------------!
    do j1=1,3
      call Get_Stat(Dist2_u(j1) , average,var,abs,min,max,N , Warning=.false.)
      d2max(j1)=TCfactor**2 * average
      call Get_Stat(Dist2_v(j1) , average,var,abs,min,max,N , Warning=.false.)
      d2max(j1)=d2max(j1) + TCfactor**2 * average
    enddo
  enddo Iteration_Loop
!------------------------------------------------------------------------------!
! Write output.                                                                !
!------------------------------------------------------------------------------!
  write (*,*) ' '
  write (*,*) '           Z O N A L   W I N D   C O M P   U       M E R I D   W I N D   C O M P   V'
  write (*,*) 'Inst      bias     scale     sigma    dsigma      bias     scale     sigma    dsigma'
  write (*,*) '------------------------------------------------------------------------------------'

  select case (Reference)
  case ('buoy')
    write (*,'(a5,8f10.6)') ' Buoy',b_u(1),a_u(1),eps_u(1),deps_u(1),  &
                                    b_v(1),a_v(1),eps_v(1),deps_v(1)
    write (*,'(a5,8f10.6)') ' Scat',b_u(2),a_u(2),eps_u(2),deps_u(2),  &
                                    b_v(2),a_v(2),eps_v(2),deps_v(2)
    write (*,'(a5,8f10.6)') ' Back',b_u(3),a_u(3),eps_u(3),deps_u(3),  &
                                    b_v(3),a_v(3),eps_v(3),deps_v(3)
  case ('scat')
    write (*,'(a5,8f10.6)') ' Buoy',b_u(3),a_u(3),eps_u(3),deps_u(3),  &
                                    b_v(3),a_v(3),eps_v(3),deps_v(3)
    write (*,'(a5,8f10.6)') ' Scat',b_u(1),a_u(1),eps_u(1),deps_u(1),  &
                                    b_v(1),a_v(1),eps_v(1),deps_v(1)
    write (*,'(a5,8f10.6)') ' Back',b_u(2),a_u(2),eps_u(2),deps_u(2),  &
                                    b_v(2),a_v(2),eps_v(2),deps_v(2)
  case ('back')
    write (*,'(a5,8f10.6)') ' Buoy',b_u(2),a_u(2),eps_u(2),deps_u(2),  &
                                    b_v(2),a_v(2),eps_v(2),deps_v(2)
    write (*,'(a5,8f10.6)') ' Scat',b_u(3),a_u(3),eps_u(3),deps_u(3),  &
                                    b_v(3),a_v(3),eps_v(3),deps_v(3)
    write (*,'(a5,8f10.6)') ' Back',b_u(1),a_u(1),eps_u(1),deps_u(1),  &
                                    b_v(1),a_v(1),eps_v(1),deps_v(1)
  case default
  end select

  end subroutine CalcTCforFile_Jur


  subroutine GetTripleCollocationWinds(Line , u,v)
!------------------------------------------------------------------------------!
! Reads the triple collocation winds from the data line and stores them in u,v.!
!------------------------------------------------------------------------------!
  implicit none

  type(DataLine_Type), intent(in)  :: Line  ! Data line struct
  real, dimension(3),intent(out)   :: u,v   ! Wind components

  select case (Reference)
  case ('buoy')
    u(1)=Line%Buoyu
    v(1)=Line%Buoyv
    u(2)=Line%Scatu
    v(2)=Line%Scatv
    u(3)=Line%Backu
    v(3)=Line%Backv
  case ('scat')
    u(1)=Line%Scatu
    v(1)=Line%Scatv
    u(2)=Line%Backu
    v(2)=Line%Backv
    u(3)=Line%Buoyu
    v(3)=Line%Buoyv
  case ('back')
    u(1)=Line%Backu
    v(1)=Line%Backv
    u(2)=Line%Buoyu
    v(2)=Line%Buoyv
    u(3)=Line%Scatu
    v(3)=Line%Scatv
  case default
    write (*,*) 'ERROR in GetTripleCollocationWinds: illegal reference ',Reference
    stop
  end select
 
  end subroutine GetTripleCollocationWinds

!==============================================================================!
! ROUTINES FOR DOUBLE COLLOCATION.                                             !
!==============================================================================!

  subroutine CalcDoubleCollocation
!------------------------------------------------------------------------------!
! Calculate double collocation.                                                !
!------------------------------------------------------------------------------!
  implicit none

  integer                        :: i

  if (MaxIter2 == 0) return

  do i=1,Set%NrOfFiles
    call CalcDCforFile(Set%File(i)  , DoubleCollocation%DC(i))
    call CalcDCforFile_New(Set%File(i))
  enddo

  end subroutine CalcDoubleCollocation


  subroutine CalcDCforFile(File,DC)
!------------------------------------------------------------------------------!
! Evaluates the double collocation error model with third order statistics.    !
!------------------------------------------------------------------------------!
  implicit none

  type(DataFile_Type), intent(inout)    :: File                         ! Data file struct
  type(Collocation_Type), intent(inout) :: DC                           ! Double collocation moments struct
!------------------------------------------------------------------------------!
! Local declarations.                                                          !
!------------------------------------------------------------------------------!
  real, dimension(2)                    :: u,v                          ! Vector of (u,v) values
  real                                  :: a_u,a_v,bx_u,bx_v,by_u,by_v  ! Calibration coefficients
  real                                  :: d2u,d2v                      ! Distance between wind components
  type(Stat_Type)                       :: Dist2_u,Dist2_v              ! Distance statistics struct
  real                                  :: d2,d2max                     ! Actual/maximum distance from calibration curve
  real                                  :: epsx_u,epsx_v,epsy_u,epsy_v  ! 2nd order error statistics (width)
  real                                  :: gamx_u,gamx_v,gamy_u,gamy_v  ! 3rd order error statistics (skewness)
  real                                  :: t2_u,t2_v                    ! Truth 2nd order variances
  real                                  :: t3_u,t3_v                    ! Truth 3rd order variances
  real                                  :: average,var,abs,min,max      ! Statistics parameters
  integer                               :: N                            ! Number of points in this file
  integer                               :: Iteration                    ! Iteration number
  integer                               :: Rejected                     ! Number of data lines rejected by sigma-test
  integer                               :: l                            ! Index
!------------------------------------------------------------------------------!
! Initialize. Start iteration.                                                 !
!------------------------------------------------------------------------------!
  bx_u=0.0
  by_u=0.0
  bx_v=0.0
  by_v=0.0

  a_u=1.0
  a_v=1.0

  d2max=TCfactor**2 * 9.0

  if (Verbosity > 1) then
    write (*,*) ' '
    write (*,*) 'Double Collocation intermediate results:'
    write (*,*) ' '

    write (*,*) '                 S C A T             B A C K'
    write (*,*) '        Error Std.Dev. 1    Error Std.Dev. 2'           &
                  //'   Scale parameter 2      Truth Std.Dev.'
    write (*,*) 'Iter         u         v         u         v'           &
                  //'         u         v         u         v  Rejected'
    write (*,*) '--------------------------------------------'           &
                  //'--------------------------------------------------'
  end if

  Iteration_Loop: do Iteration=1,MaxIter2
!------------------------------------------------------------------------------!
!   Reset moments. Start loop over data points.                                !
!------------------------------------------------------------------------------!
    call Reset_Moments(DC%u)
    call Reset_Moments(DC%v)

    call Init_Stat(Dist2_u)
    call Init_Stat(Dist2_v)

    Rejected=0

    Data_Loop: do l=1,File%NrOfDataLines
      u(1)=File%Line(l)%Scatu
      u(2)=File%Line(l)%Backu
      v(1)=File%Line(l)%Scatv
      v(2)=File%Line(l)%Backv

!!    u(1)=File%Line(l)%Backu
!!    u(2)=File%Line(l)%Scatu
!!    v(1)=File%Line(l)%Backv
!!    v(2)=File%Line(l)%Scatv
!------------------------------------------------------------------------------!
!     Sigma test.                                                              !
!------------------------------------------------------------------------------!
      d2u=(bx_u + u(1) - by_u - a_u*u(2))**2
      d2v=(bx_v + v(1) - by_v - a_v*v(2))**2

      d2=d2u + d2v
      if (d2 >= d2max) then
        Rejected=Rejected + 1
        cycle Data_Loop
      end if
!------------------------------------------------------------------------------!
!     Update first and second moments; update distances                        !
!------------------------------------------------------------------------------!
      call Update_Moments(DC%u , u)
      call Update_Moments(DC%v , v)

      call Update_Stat(Dist2_u,d2u)
      call Update_Stat(Dist2_v,d2v)
    enddo Data_Loop
!------------------------------------------------------------------------------!
!   Calculate error model parameters.                                          !
!------------------------------------------------------------------------------!
    call DoDC(DC%u , t3_u,t2_u , gamx_u,gamy_u , epsx_u,epsy_u , a_u,bx_u,by_u)
    call DoDC(DC%v , t3_v,t2_v , gamx_v,gamy_v , epsx_v,epsy_v , a_v,bx_v,by_v)

!!  call PrintDCMoments

    if (Verbosity > 1) write (*,'(i5,8f10.5,i10)') Iteration,epsx_u,epsx_v,epsy_u,epsy_v, &
                                                   a_u,a_v,sqrt(t2_u),sqrt(t2_v),Rejected
!------------------------------------------------------------------------------!
!   Prepare next iteration.                                                    !
!------------------------------------------------------------------------------!
    call Get_Stat(Dist2_u , average,var,abs,min,max,N , Warning=.false.)
    d2max=TCfactor**2 * average
    call Get_Stat(Dist2_v , average,var,abs,min,max,N , Warning=.false.)
    d2max=d2max + TCfactor**2 * average
  enddo Iteration_Loop

  end subroutine CalcDCforFile


  subroutine CalcDCforFile_New(File)
!------------------------------------------------------------------------------!
! Calculates for a single file the moments for double collocation using only   !
! those points that are not too far from the calibration line.                 !
! This is  Ad stoffelen's formulation with optimized statistics.               !
!------------------------------------------------------------------------------!
  implicit none

  type(DataFile_Type), intent(inout)    :: File                     ! Data file struct
!------------------------------------------------------------------------------!
! Local declarations.                                                          !
!------------------------------------------------------------------------------!
  real                                  :: Mx_u,Mx_v,My_u,My_v      ! First order moments
  real                                  :: Dx_u,Dx_v,Dy_u,Dy_v      ! Second order moments
  real                                  :: Mxy_u,Mxy_v              ! Second order moments
  real                                  :: Dxx_u,Dxx_v,Dyy_u,Dyy_v  ! Third order moments
  real                                  :: Mxxy_u,Mxxy_v            ! Third order moments
  real                                  :: Mxyy_u,Mxyy_v            ! Third order moments

  real, dimension(2)                    :: uin,vin                  ! Vector of input (u,v) values
  real, dimension(2)                    :: u,v                      ! Vector of calibrated (u,v) values
  real, dimension(2)                    :: Bias_u,Bias_v            ! Calibration coefficients
  real, dimension(2)                    :: Scale_u,Scale_v          ! Calibration coefficients
  real, dimension(2)                    :: BiasCorrection_u
  real, dimension(2)                    :: BiasCorrection_v
  real, dimension(2)                    :: ScaleCorrection_u
  real, dimension(2)                    :: ScaleCorrection_v
  real                                  :: d2u,d2v                  ! Distance between wind components
  type(Stat_Type)                       :: Dist2_u,Dist2_v          ! Distance statistics struct
  real                                  :: d2,d2max                 ! Actual/maximum distance from calibration curve
  real, dimension(2)                    :: eps_u,eps_v              ! Triple collocation error standard deviation
  real, dimension(2)                    :: deps_u,deps_v            ! Precision in triple collocation error standard deviation
  real, dimension(2)                    :: gam_u,gam_v
  real                                  :: average,var,abs,min,max  ! Statistics parameters
  real                                  :: RN
  integer                               :: N                        ! Number of points in this file
  integer                               :: Iteration                ! Iteration number
  integer                               :: Rejected                 ! Number of data lines rejected by sigma-test
  integer                               :: j1,j2,l                  ! Indices
!------------------------------------------------------------------------------!
! Initialize. Start iteration.                                                 !
!------------------------------------------------------------------------------!
  do j1=1,2
    Bias_u(j1)=0.0
    Bias_v(j1)=0.0

    BiasCorrection_u(j1)=0.0
    BiasCorrection_v(j1)=0.0

    Scale_u(j1)=1.0
    Scale_v(j1)=1.0

    ScaleCorrection_u(j1)=1.0
    ScaleCorrection_v(j1)=1.0
  enddo

  d2max=TCfactor**2 * 9.0

  if (Verbosity > 1) then
    write (*,*) ' '
    write (*,*) 'Double Collocation intermediate results:'
    write (*,*) 'Method: calibrated data'
    write (*,*) ' '

    write (*,*) '                 S C A T             B A C K'
    write (*,*) '        Error Std.Dev. 1    Error Std.Dev. 2'           &
                  //'   Scale parameter 2'
    write (*,*) 'Iter         u         v         u         v'           &
                  //'         u         v  Rejected'
    write (*,*) '--------------------------------------------'           &
                  //'------------------------------'
  end if

  Iteration_Loop: do Iteration=1,MaxIter2
!------------------------------------------------------------------------------!
!   Reset moments. Start loop over data points.                                !
!------------------------------------------------------------------------------!
    Mx_u=0.0
    Mx_v=0.0
    My_u=0.0
    My_v=0.0

    Mxy_u=0.0
    Mxy_v=0.0
    Dx_u=0.0
    Dx_v=0.0
    Dy_u=0.0
    Dy_v=0.0

    Mxxy_u=0.0
    Mxxy_v=0.0
    Mxyy_u=0.0
    Mxyy_v=0.0
    Dxx_u=0.0
    Dxx_v=0.0
    Dyy_u=0.0
    Dyy_v=0.0
    
    call Init_Stat(Dist2_u)
    call Init_Stat(Dist2_v)

    Rejected=0
    N=0

    Data_Loop: do l=1,File%NrOfDataLines
      uin(1)=File%Line(l)%Scatu
      uin(2)=File%Line(l)%Backu
      vin(1)=File%Line(l)%Scatv
      vin(2)=File%Line(l)%Backv
!------------------------------------------------------------------------------!
!     Calibration. BiasCorrection remains zero if -nobias option is set.       !
!------------------------------------------------------------------------------!
      do j1=1,2
        u(j1)=(uin(j1) - BiasCorrection_u(j1))*ScaleCorrection_u(j1)
        v(j1)=(vin(j1) - BiasCorrection_v(j1))*ScaleCorrection_v(j1)
      enddo
!------------------------------------------------------------------------------!
!     Sigma test.                                                              !
!------------------------------------------------------------------------------!
      d2u=(u(1) - u(2))**2
      d2v=(v(1) - v(2))**2

      d2=d2u + d2v
      if (d2 >= d2max) then
        Rejected=Rejected + 1
        cycle Data_Loop
      end if
!------------------------------------------------------------------------------!
!     Update optimized second moments and distances                            !
!------------------------------------------------------------------------------!
      N=N + 1
      RN=1.0/real(N)

      Mx_u=Mx_u + RN*(u(1) - Mx_u)
      Mx_v=Mx_v + RN*(v(1) - Mx_v)
      My_u=My_u + RN*(u(2) - My_u)
      My_v=My_v + RN*(v(2) - My_v)

      Mxy_u=Mxy_u + RN*(u(1)*u(2) - Mxy_u)
      Mxy_v=Mxy_v + RN*(v(1)*v(2) - Mxy_v)
      Dx_u=Dx_u + RN*(u(1)*(u(1) - u(2)) - Dx_u)
      Dx_v=Dx_v + RN*(v(1)*(v(1) - v(2)) - Dx_v)
      Dy_u=Dy_u + RN*(u(2)*(u(1) - u(2)) - Dy_u)
      Dy_v=Dy_v + RN*(v(2)*(v(1) - v(2)) - Dy_v)

      Mxxy_u=Mxxy_u + RN*(u(1)*u(1)*u(2) - Mxxy_u)
      Mxxy_v=Mxxy_v + RN*(v(1)*v(1)*v(2) - Mxxy_v)
      Mxyy_u=Mxyy_u + RN*(u(1)*u(2)*u(2) - Mxyy_u)
      Mxyy_v=Mxyy_v + RN*(v(1)*v(2)*v(2) - Mxyy_v)
      Dxx_u=Dxx_u + RN*(u(1)*u(1)*(u(1) - u(2)) - Dxx_u)
      Dxx_v=Dxx_v + RN*(v(1)*v(1)*(v(1) - v(2)) - Dxx_v)
      Dyy_u=Dyy_u + RN*(u(2)*u(2)*(u(1) - u(2)) - Dyy_u)
      Dyy_v=Dyy_v + RN*(v(2)*v(2)*(v(1) - v(2)) - Dyy_v)

      call Update_Stat(Dist2_u,d2u)
      call Update_Stat(Dist2_v,d2v)
    enddo Data_Loop
!------------------------------------------------------------------------------!
!   Calculate double collocation errors.                                       !
!------------------------------------------------------------------------------!
    eps_u(1)=sqrt(Dx_u)
    eps_u(2)=sqrt(Dy_v)
    eps_v(1)=sqrt(Dx_v)
    eps_v(2)=sqrt(Dy_v)

    gam_u(1)=Dxx_u
    gam_u(2)=Dyy_u
    gam_v(1)=Dxx_v
    gam_v(2)=Dyy_v

    if (Verbosity > 1) write (*,'(i5,6f10.5,i10)') Iteration,eps_u(1),eps_v(1),eps_u(2),eps_v(2), &
                                   Scale_u(2),Scale_v(2),Rejected
!------------------------------------------------------------------------------!
!   Calculate precision in triple collocation error parameters.                !
!------------------------------------------------------------------------------!
    if (Iteration == MaxIter2) then
      RN=1.0/sqrt(real(N))

      do j1=1,2
        j2=mod(j1,2) + 1
        deps_u(j1)=RN*eps_u(j1)*sqrt(2.0*eps_u(j1)**2 + eps_u(j2)**2)
        deps_v(j1)=RN*eps_v(j1)*sqrt(2.0*eps_v(j1)**2 + eps_v(j2)**2)
      enddo

      if (Verbosity > 1) write (*,'(a5,4f10.5)') ' Prec',deps_u(1),deps_v(1),deps_u(2),deps_v(2)
    end if
!------------------------------------------------------------------------------!
!   Prepare next iteration.                                                    !
!------------------------------------------------------------------------------!
    if (DoBiasCorrection) then
      Bias_u(1)=Mx_u
      Bias_u(2)=My_u
      Bias_v(1)=Mx_v
      Bias_v(2)=My_v
    end if

    Scale_u(1)=1.0
    Scale_u(2)=Dxx_u/Dyy_u
    Scale_v(1)=1.0
    Scale_v(2)=Dxx_v/Dyy_v

    do j1=1,2
      if (DoBiasCorrection) then
        BiasCorrection_u(j1)=BiasCorrection_u(j1) + Bias_u(j1)
        BiasCorrection_v(j1)=BiasCorrection_v(j1) + Bias_v(j1)
      end if

      ScaleCorrection_u(j1)=ScaleCorrection_u(j1)/Scale_u(j1)
      ScaleCorrection_v(j1)=ScaleCorrection_v(j1)/Scale_v(j1)
    enddo

    call Get_Stat(Dist2_u , average,var,abs,min,max,N , Warning=.false.)
    d2max=TCfactor**2 * average
    call Get_Stat(Dist2_v , average,var,abs,min,max,N , Warning=.false.)
    d2max=d2max + TCfactor**2 * average
  enddo Iteration_Loop

  end subroutine CalcDCforFile_New


  subroutine DoDC(C , t3,t2 , gamx,gamy , epsx,epsy , a,bx,by)
!------------------------------------------------------------------------------!
! Performs the double collocation error model calculation.                     !
!------------------------------------------------------------------------------!
  implicit none

  type (Moments_Type), intent(in)  :: C
  real, intent(out)                :: t3,t2
  real, intent(out)                :: gamx,gamy
  real, intent(out)                :: epsx,epsy
  real, intent(out)                :: a,bx,by

  real                             :: Cxx,Cxy,Cyy
  real                             :: Dxxx,Dxxy,Dxyy,Dyyy
!------------------------------------------------------------------------------!
! Set biases. Set covariances and 3rd order quantities.                        !
!------------------------------------------------------------------------------!
  bx=C%M1(1)
  by=C%M1(2)

  Cxx=C%M2(1,1) - bx*bx
  Cxy=C%M2(1,2) - bx*by
  Cyy=C%M2(2,2) - by*by

  Dxxx=C%M3(1,1,1) - 3.0*Cxx*bx - bx**3
  Dxxy=C%M3(1,1,2) - Cxx*by - 2.0*bx*Cxy - bx**2*by
  Dxyy=C%M3(1,2,2) - 2.0*Cxy*by - bx*Cyy - bx*by**2
  Dyyy=C%M3(2,2,2) - 3.0*Cyy*by - by**3

!!write (*,*) ' '
!!write (*,'(a,3f12.6)') 'DoDC: Cxx,Cxy,Cyy        =',Cxx,Cxy,Cyy
!!write (*,'(a,4f12.6)') 'DoDC: Dxxx,Dxxy,Dxyy,Dyyy=',Dxxx,Dxxy,Dxyy,Dyyy
!------------------------------------------------------------------------------!
! Calculate error model parameters.                                            !
!------------------------------------------------------------------------------!
  t3=Dxxy**2/Dxyy

  a=Dxyy/Dxxy

  gamx=Dxxx - t3
  gamy=Dyyy - a**3*t3

  t2=Cxy/a

  epsx=sqrt(Cxx - t2)
  epsy=sqrt(Cyy - a**2*t2)

  end subroutine DoDC

!==============================================================================!
! ROUTINES FOR OUTPUT ON FILE.                                                 !
!==============================================================================!

  subroutine SaveSet
!------------------------------------------------------------------------------!
! Writes the data set to file.                                                 !
!------------------------------------------------------------------------------!
  implicit none

  character(len=256)   :: FileName
  integer              :: file,line

  write (*,*) ' '
  write (*,*) 'Results of routine SaveSet:'
  write (*,*) '  File  Action'
  write (*,*) '  --------------------------------'

  do file=1,Set%NrOfFiles
    FileName=trim(Set%File(file)%Name) // "." // trim(Extension)
    open (saveLUN, file=FileName, status='new')

    write (*,'(i7,2x,a)') file,'Saved to '//trim(FileName)

    do line=1,Set%File(file)%NrOfDataLines
      write (saveLUN,'(i5,2(2f8.2,i10,i5,2f6.1),2f6.1)')       &
                        Set%File(file)%Line(line)%BuoyNumber,  &
                        Set%File(file)%Line(line)%BuoyLat,     &
                        Set%File(file)%Line(line)%BuoyLon,     &
                        Set%File(file)%Line(line)%BuoyDate,    &
                        Set%File(file)%Line(line)%BuoyTime,    &
                        Set%File(file)%Line(line)%BuoySpeed,   &
                        Set%File(file)%Line(line)%BuoyDir,     &
                        Set%File(file)%Line(line)%ScatLat,     &
                        Set%File(file)%Line(line)%ScatLon,     &
                        Set%File(file)%Line(line)%ScatDate,    &
                        Set%File(file)%Line(line)%ScatTime,    &
                        Set%File(file)%Line(line)%ScatSpeed,   &
                        Set%File(file)%Line(line)%ScatDir,     &
                        Set%File(file)%Line(line)%BackSpeed,   &
                        Set%File(file)%Line(line)%BackDir
    enddo
  enddo

  end subroutine SaveSet

!==============================================================================!
! ROUTINES FOR OUTPUT ON SCREEN.                                               !
!==============================================================================!

  subroutine DumpFile(File)
!------------------------------------------------------------------------------!
! Writes the contents of a File to screen.                                     !
!------------------------------------------------------------------------------!
  implicit none

  type(DataFile_Type), intent(in)  :: File
  integer                          :: l

  write (*,*) 'Contents of file ',trim(File%Name)
  write (*,*) '         BUOY                                               '// &
              'SCATTEROMETER                         BACKGROUND'
  write (*,*) '  Line   Buoy     Lat     Lon     Date Time Speed Direc     '// &
              'Lat     Lon     Date Time Speed Direc Speed Direc  Common'
  write (*,*) '------------------------------------------------------------'// &
              '---------------------------------------------------------'
  do l=1,File%NrOfDataLines
    write (*,'(2i7,2(2f8.2,i9,i5,2f6.1),2f6.1,l8)') l ,File%Line(l)%BuoyNumber, &
      File%Line(l)%BuoyLat   , File%Line(l)%BuoyLon   ,File%Line(l)%BuoyDate,   &
      File%Line(l)%BuoyTime  , File%Line(l)%BuoySpeed ,File%Line(l)%BuoyDir,    &
      File%Line(l)%ScatLat   , File%Line(l)%ScatLon   ,File%Line(l)%ScatDate,   &
      File%Line(l)%ScatTime  , File%Line(l)%ScatSpeed ,File%Line(l)%ScatDir,    &
      File%Line(l)%BackSpeed , File%Line(l)%BackDir   ,File%Line(l)%Common
  enddo

  end subroutine DumpFile


  subroutine PrintCollocationStatistics
!------------------------------------------------------------------------------!
! Prints the results for the collocation statistics.                           !
!------------------------------------------------------------------------------!
  implicit none

  character(len=72) :: Header1,Header2

  Header1='Results for Scatterometer minus buoy'
  Header2='Results for Scatterometer minus background'

  write (*,*) ' '
  write (*,*) ' '
  write (*,*) '=================================='
  write (*,*) 'STATISTICS FOR COMMON COLLOCATIONS'
  write (*,*) '=================================='
  call PrintCSStruct(CommonStatistics_ScatBuoy,Header1)
  call PrintCSStruct(CommonStatistics_ScatBack,Header2)

  write (*,*) ' '
  write (*,*) ' '
  write (*,*) '======================================'
  write (*,*) 'STATISTICS FOR INDIVIDUAL COLLOCATIONS'
  write (*,*) '======================================'
  call PrintCSStruct(IndividualStatistics_ScatBuoy,Header1)
  call PrintCSStruct(IndividualStatistics_ScatBack,Header2)

  write (*,*) ' '
  write (*,*) ' '
  write (*,*) '======================================'
  write (*,*) 'STATISTICS FOR NON-COMMON COLLOCATIONS'
  write (*,*) '======================================'
  call PrintCSStruct(NonCommonStatistics_ScatBuoy,Header1)
  call PrintCSStruct(NonCommonStatistics_ScatBack,Header2)

  end subroutine PrintCollocationStatistics


  subroutine PrintCSStruct(CollStat,Header)
!------------------------------------------------------------------------------!
! Prints the results for a single collocation statistics struct.               !
!------------------------------------------------------------------------------!
  implicit none

  type(CollocationStatistics_Type), intent(in) :: CollStat
  character(len=72), intent(in)                :: Header
  real                                         :: Bias,Var,StdDev,Abs,Min,Max
  integer                                      :: N
  integer                                      :: i

  write (*,*) ' '
  write (*,*) trim(Header)
  write (*,*) 'Results for difference in wind speed:'
  write (*,*) 'File    Bias     Var  StdDev     Abs     Min     Max     N  Name'
  write (*,*) '----------------------------------------------------------------------------------'
  do i=1,Set%NrOfFiles
    call Get_Stat(CollStat%S(i) , Bias,Var,Abs,Min,Max,N , Warning=.false.)
    StdDev=0.0
    if (Var > 0.0) StdDev=sqrt(Var)
    write (*,'(i5,4f8.3,2f8.2,i6,2x,a)') i,Bias,Var,StdDev,Abs,Min,Max,N,trim(Set%File(i)%Name)
  enddo

  write (*,*) ' '
  write (*,*) trim(Header)
  write (*,*) 'Results for difference in wind direction:'
  write (*,*) 'File    Bias     Var  StdDev     Abs     Min     Max     N  Name'
  write (*,*) '----------------------------------------------------------------------------------'
  do i=1,Set%NrOfFiles
    call Get_Stat(CollStat%D(i) , Bias,Var,Abs,Min,Max,N , Warning=.false.)
    StdDev=0.0
    if (Var > 0.0) StdDev=sqrt(Var)
    write (*,'(i5,4f8.3,2f8.2,i6,2x,a)') i,Bias,Var,StdDev,Abs,Min,Max,N,trim(Set%File(i)%Name)
  enddo

  write (*,*) ' '
  write (*,*) trim(Header)
  write (*,*) 'Results for difference in zonal wind component u:'
  write (*,*) 'File    Bias     Var  StdDev     Abs     Min     Max     N  Name'
  write (*,*) '----------------------------------------------------------------------------------'
  do i=1,Set%NrOfFiles
    call Get_Stat(CollStat%u(i) , Bias,Var,Abs,Min,Max,N , Warning=.false.)
    StdDev=0.0
    if (Var > 0.0) StdDev=sqrt(Var)
    write (*,'(i5,4f8.3,2f8.2,i6,2x,a)') i,Bias,Var,StdDev,Abs,Min,Max,N,trim(Set%File(i)%Name)
  enddo

  write (*,*) ' '
  write (*,*) trim(Header)
  write (*,*) 'Results for difference in meridional wind component v:'
  write (*,*) 'File    Bias     Var  StdDev     Abs     Min     Max     N  Name'
  write (*,*) '----------------------------------------------------------------------------------'
  do i=1,Set%NrOfFiles
    call Get_Stat(CollStat%v(i) , Bias,Var,Abs,Min,Max,N , Warning=.false.)
    StdDev=0.0
    if (Var > 0.0) StdDev=sqrt(Var)
    write (*,'(i5,4f8.3,2f8.2,i6,2x,a)') i,Bias,Var,StdDev,Abs,Min,Max,N,trim(Set%File(i)%Name)
  enddo

  end subroutine PrintCSStruct


  subroutine PrintTCMoments
!------------------------------------------------------------------------------!
! Prints the first and second moments used for triple collocation.             !
!------------------------------------------------------------------------------!
  implicit none

  integer  :: i

  write (*,*) ' '
  write (*,*) 'First moments used for triple collocation'
  write (*,*) '                    zonal wind component u         meridional wind component v'
  write (*,*) '  File          M1          M2          M3          M1          M2          M3       N'
  write (*,*) '  ------------------------------------------------------------------------------------'

  do i=1,TripleCollocation%N
    write (*,'(i7,6f12.6,i8)') i,TripleCollocation%TC(i)%u%M1(1),TripleCollocation%TC(i)%u%M1(2),  &
                                 TripleCollocation%TC(i)%u%M1(3),TripleCollocation%TC(i)%v%M1(1),  &
                                 TripleCollocation%TC(i)%v%M1(2),TripleCollocation%TC(i)%v%M1(3),  &
                                 TripleCollocation%TC(i)%v%NrOfPoints
  enddo

  write (*,*) ' '
  write (*,*) 'Second moments of zonal wind component u used for triple collocation'
  write (*,*) '  File         M11         M22         M33         M12         M23         M31       N'
  write (*,*) '  ------------------------------------------------------------------------------------'
  do i=1,TripleCollocation%N
    write (*,'(i7,6f12.6,i8)') i,TripleCollocation%TC(i)%u%M2(1,1),TripleCollocation%TC(i)%u%M2(2,2),  &
                                 TripleCollocation%TC(i)%u%M2(3,3),TripleCollocation%TC(i)%u%M2(1,2),  &
                                 TripleCollocation%TC(i)%u%M2(2,3),TripleCollocation%TC(i)%u%M2(3,1),  &
                                 TripleCollocation%TC(i)%u%NrOfPoints
  enddo

  write (*,*) ' '
  write (*,*) 'Second moments of meridional wind component v used for triple collocation'
  write (*,*) '  File         M11         M22         M33         M12         M23         M31       N'
  write (*,*) '  ------------------------------------------------------------------------------------'
  do i=1,TripleCollocation%N
    write (*,'(i7,6f12.6,i8)') i,TripleCollocation%TC(i)%v%M2(1,1),TripleCollocation%TC(i)%v%M2(2,2),  &
                                 TripleCollocation%TC(i)%v%M2(3,3),TripleCollocation%TC(i)%v%M2(1,2),  &
                                 TripleCollocation%TC(i)%v%M2(2,3),TripleCollocation%TC(i)%v%M2(3,1),  &
                                 TripleCollocation%TC(i)%v%NrOfPoints
  enddo

  end subroutine PrintTCMoments


  subroutine PrintDCMoments
!------------------------------------------------------------------------------!
! Prints the first, second, and third moments used for double collocation.     !
!------------------------------------------------------------------------------!
  implicit none

  integer  :: i

  write (*,*) ' '
  write (*,*) 'Moments of u used for double collocation:'
  write (*,*) '  File          Mx          My         Mxx         Mxy         Myy' //  &
                    '        Mxxx        Mxxy        Mxyy        Myyy       N'
  write (*,*) '  ----------------------------------------------------------------' //  &
                    '--------------------------------------------------------'
  do i=1,DoubleCollocation%N
    write (*,'(i7,9f12.6,i8)') i,DoubleCollocation%DC(i)%u%M1(1)    ,DoubleCollocation%DC(i)%u%M1(2),      &
                                 DoubleCollocation%DC(i)%u%M2(1,1)  ,DoubleCollocation%DC(i)%u%M2(1,2),    &
                                 DoubleCollocation%DC(i)%u%M2(2,2),                                        &
                                 DoubleCollocation%DC(i)%u%M3(1,1,1),DoubleCollocation%DC(i)%u%M3(1,1,2),  &
                                 DoubleCollocation%DC(i)%u%M3(1,2,2),DoubleCollocation%DC(i)%u%M3(2,2,2),  &
                                 DoubleCollocation%DC(i)%u%NrOfPoints
  enddo

  write (*,*) ' '
  write (*,*) 'Moments of v used for double collocation:'
  write (*,*) '  File          Mx          My         Mxx         Mxy         Myy' //  &
                    '        Mxxx        Mxxy        Mxyy        Myyy       N'
  write (*,*) '  ----------------------------------------------------------------' //  &
                    '--------------------------------------------------------'
  do i=1,DoubleCollocation%N
    write (*,'(i7,9f12.6,i8)') i,DoubleCollocation%DC(i)%v%M1(1)    ,DoubleCollocation%DC(i)%v%M1(2),      &
                                 DoubleCollocation%DC(i)%v%M2(1,1)  ,DoubleCollocation%DC(i)%v%M2(1,2),    &
                                 DoubleCollocation%DC(i)%v%M2(2,2),                                        &
                                 DoubleCollocation%DC(i)%v%M3(1,1,1),DoubleCollocation%DC(i)%v%M3(1,1,2),  &
                                 DoubleCollocation%DC(i)%v%M3(1,2,2),DoubleCollocation%DC(i)%v%M3(2,2,2),  &
                                 DoubleCollocation%DC(i)%u%NrOfPoints 
  enddo

  end subroutine PrintDCMoments
  end module CollocationStatistics_Module


  program CollocationStatistics
!------------------------------------------------------------------------------!
! Main program.                                                                !
!------------------------------------------------------------------------------!
  use CollocationStatistics_Module, only: ReadArguments,CheckArguments,PrintArguments
  use CollocationStatistics_Module, only: InitSet,ReadSet,ExitSet
  use CollocationStatistics_Module, only: InitCollocationStatistics,ExitCollocationStatistics
  use CollocationStatistics_Module, only: FindCommonLines
  use CollocationStatistics_Module, only: CalcCollocationStatistics,PrintCollocationStatistics
  use CollocationStatistics_Module, only: InitTripleCollocation,ExitTripleCollocation
  use CollocationStatistics_Module, only: CalcTripleCollocation
  use CollocationStatistics_Module, only: InitDoubleCollocation,ExitDoubleCollocation
  use CollocationStatistics_Module, only: CalcDoubleCollocation
  use CollocationStatistics_Module, only: NrOfFiles
  use CollocationStatistics_Module, only: Reference

  implicit none
  character(len=4) :: Ref
!------------------------------------------------------------------------------!
! Read command line arguments. Read data.                                      !
!------------------------------------------------------------------------------!
  call ReadArguments
  call CheckArguments
  call PrintArguments

  call InitSet
  call ReadSet
!------------------------------------------------------------------------------!
! Statistics of differences between buoy and scatterometer wind.               !
!------------------------------------------------------------------------------!
  if (NrOfFiles > 1) call FindCommonLines

  call InitCollocationStatistics
  call CalcCollocationStatistics
  call PrintCollocationStatistics
  call ExitCollocationStatistics
!------------------------------------------------------------------------------!
! Triple collocation.                                                          !
!------------------------------------------------------------------------------!
  Ref=Reference                 ! Make a copy because Reference may be changed
  write (*,*) 'Reference for TC is ',Ref
  select case (Ref)
  case ('all ' , ' all')
    write (*,*) ' '
    write (*,*) ' '
    write (*,*) '=========================='
    write (*,*) 'TRIPLE COLLOCATION RESULTS'
    write (*,*) '=========================='

    Reference='buoy'
    call InitTripleCollocation
    call CalcTripleCollocation
    call ExitTripleCollocation

    Reference='scat'
    call InitTripleCollocation
    call CalcTripleCollocation
    call ExitTripleCollocation

    Reference='back'
    call InitTripleCollocation
    call CalcTripleCollocation
    call ExitTripleCollocation
  case ('buoy' , 'scat' , 'back')
    write (*,*) ' '
    write (*,*) ' '
    write (*,*) '=========================='
    write (*,*) 'TRIPLE COLLOCATION RESULTS'
    write (*,*) '=========================='

    call InitTripleCollocation
    call CalcTripleCollocation
    call ExitTripleCollocation
  case ('none')
  case default
  end select
!------------------------------------------------------------------------------!
! Double collocation with third-order statistics.                              !
!------------------------------------------------------------------------------!
  call InitDoubleCollocation
  call CalcDoubleCollocation
  call ExitDoubleCollocation
!------------------------------------------------------------------------------!
! Finalise.                                                                    !
!------------------------------------------------------------------------------!
  call ExitSet

  end program CollocationStatistics


